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In the FORTRAN notation in a write statement 14 lines from the bottom of the 
page, H should be changed to AH. The corrected notation would be as 
follows: 


WRITE(6,126)DTTI,DTTF,DTWI,DTWF,RMIN,AH 

Insert the following statements after FORTRAN statement 402 (before 8th line 
from bottom of page): 

C TEST FOR ZERO VALUED DENOMINATOR 

CKDEN=RHO l*COS(TH 1) -RH02*C0S(TH2) 

IF(CKDEN.GE.1.0E- 10)410, 405 
410 CONTINUE 

Insert the following statements after FORTRAN statement 602 (after 7th line 
from top of page): 

C TEST FOR ZERO VALUED DENOMINATOR 

CKDEN=RHO l*COS(TH 1) -RH02*C0S(TH2) 
IF(CKDEN.GE.1.0E-10)614,604 
614 CONTINUE 

The absolute value of ij and i should be used in Subroutine Six which evalu- 
ates the transformation matrix. Therefore, the second statement on page 72 
should be deleted and replaced by the following two statements: 


RAH=AH 

CALL SIX(All,A12,A21,A22,A31,A32,OMI,OMEGAI,RAII) 

Also, the present one-line FORTRAN statement 900 should be deleted and 
replaced by the following two statements: 

900 RBI=BI 

CALL SIX(Bll,Bl2,B21,B22,B31,B32,OM, OMEGA, RBI) 
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Page 78: Add the following statement between lines 1 and 2: 


U3=ABS(U3) 

Page 83 (table HI): Under the heading DV1Z (12th column), delete the minus signs 
before the first six numbers and insert minus signs before the last nine 
numbers; under the heading DV2Z (15th column), insert minus signs 
before all fifteen numbers. 
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By Robert L. Collins and Sylvia A. Wallace 
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SUMMARY 

A detailed derivation of exact equations and an associated computer program with 
which the basic parameters involved in transfer or rendezvous between two arbitrary 
elliptic orbits may be calculated are presented. The computations are exact in the sense 
that the Kepler solutions for bodies orbiting in a point- gravity field are used. This 
method deviates from the "classical" approach based on the theorem of Lambert, inas- 
much as it uses the true anomaly and the Kepler equations for iterating to the desired 
rendezvous transfer time. The method has a unique feature in that definite boundaries, 
dependent on the problem input, are used which limit the range of the true anomaly, and 
thereby reduce the search effort required in the iteration procedure. The program pro- 
vides for a solution to problems where the transfer angle is less than 360°. Examples 
are given for three particular uses of the program: (1) interplanetary transfer between 
massless planets, (2) near-planet orbit rendezvous, and (3) orbital transfer. 

INTRODUCTION 

It seems appropriate to provide a useful technique for the computation of velocity 
increments and other important parameters involved in the problems of orbital rendezvous 
and transfer by use of the solution to the exact equations of motion. Various programs 
and techniques exist at this time; however, they are either unpublished or are not of a 
sufficiently general nature to be used in the variety of orbital problems one might desire. 
(For instance, see refs. 1 to 4.) It is desirable, therefore, to have a simple and yet gen- 
eral computational method which will solve the problems of: (a) interplanetary transfer 
between massless planets, (b) planet orbit rendezvous, and (c) orbital transfer. The pur- 
pose of this paper is to give the description of the analysis required for the solution of 
the rendezvous problem from the exact Keplerian relations. The fundamental problem is 
that of determining the velocity increment required to rendezvous from some initial inter- 
ceptor orbit to some final target orbit as a function of transfer time. Also determined are 
eccentricity, semi-major axis, initial and final anomalies, and other parameters associ- 
ated with the transfer orbit. Although the problem, as stated, is a rendezvous problem, it 



is also possible to interpret the results for use in studying the orbital transfer problem. 

It is assumed that the Keplerian orbital quantities are known in advance for the target 
vehicle and that the initial interceptor orbit is known either from its Keplerian orbital 
elements or from relative coordinate data. The coordinates chosen are referenced to the 
target orbit, and another axis transformation will be necessary if the user desires the 
results and input referenced to some other axis system (such as the ecliptic). Figures 1, 
2, and 3 show the coordinates and position symbols used to describe the orbits and fig- 
ures 4 and 5 show the (input) information needed on the position of the orbits. 

A solution to the rendezvous problem may be obtained by specifying the transfer 
time for the interceptor to travel from its initial orbit to its final (target) orbit which, 
along with the initial conditions of the problem, gives the initial and final positions in 
space through which the interceptor must pass. The transfer orbit must then be found 
which passes through the known initial and final positions. An iteration is required for 
this calculation inasmuch as the transfer orbit which will yield the proper transfer time 
is not yet known. The procedure begins by choosing some orbit which passes through the 
initial and final position vectors and then the corresponding transfer time is computed 
and compared with the real (desired) transfer time. If the computed transfer time is not 
the same as the desired time, another orbit must be chosen and the time again computed 
and compared. This process is repeated until the conic section which provides the true 
rendezvous transfer orbit is found. 

There are (at least) two methods which have been used for this iteration. The solu- 
tion by use of Lambert's theorem was used by Battin (ref. 1) and Breakwell (ref. 2). How- 
ever, iteration of the true anomaly, as was considered by Lascody (ref. 3), is more read- 
ily visualized since it does not require such artificial devices as "flattening" the transfer 
orbit used in geometrically describing the transfer problem deduced from the theorem of 
Lambert (ref. 1). It is also possible to find definite regions within which the true anomaly 
must be chosen and therefore shorten computational (convergence) time to the extent that 
a direct search routine may be used in the iteration; this routine does not involve the 
derivatives necessary in the method of iteration as presented in reference 3. 

In programing this problem, an effort has been made to provide a sufficiently gen- 
eral program with a minimum amount of difficulty in reading input and printing output 
information for the particular problem. For instance, problem input information may be 
provided in two fundamentally different ways: (1) as Keplerian data (eccentricity, semi- 
major axis, longitude of ascending node, etc.) and (2) as relative coordinate input refer- 
enced to rectangular coordinate axes fixed with the target vehicle. 
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SYMBOLS 


Unless otherwise noted all quantities are nondimensional. The nondimensional 
forms are derived as follows: 


Nondimensional length = 


Nondimensional velocity = 


Dimensional length 


Semi- major axis of target orbit, a T 


Dimensional velocity 


Nondimensional angular rate = (Dimensional angular rate) x 


Mean circular velocity of target, 

Target orbital period, P T 


2ir 


a ij 


b ii 


E 

F 

F 1’ F 2 

h 

H 

i 

A A A 
A A. A 

I,J,K 


Nondimensional time = 


Dimensional time 

Target orbital period, Prp 


semi-major axis of orbit 

element of transformation matrix, x,y,z to X,Y,Z 

element of transformation matrix, x^y^z' to X,Y,Z 

chord joining p^ to P 2 

eccentricity of orbit 

eccentric anomaly of vehicle in orbit 

functional relationship 

equivalent eccentric anomalies for hyperbolic orbit 
increment for iteration 
angular momentum 

inclination of orbital plane to target orbit plane 
unit vectors along x,y,z; x",y",z M axes as indicated 
unit vectors along x'jy'jZ' axes 
unit vectors along X,Y,Z axes 
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M 

P 


mean anomaly 
semi-latus rectum 


FJ 

P T target orbital period, 27m ~ jj ~> seconds 

r radial distance to target 

r f (X),r f (Y),r f (Z) velocity components after final impulse 

t time 

t time at observation of input 

time of initial impulse 
time of final impulse 
wait time before initial impulse, t ^ - t Q 
interceptor transfer time, t f - tj 


*i 


At. 


w 


computed times from periapse when the interceptor is at p^ and p 2 in 
transfer orbit 

computed transfer time, Tg - for comparison 
speed 

mean circular speed of target orbit, if^-, ft/sec 


At 

T 1> T 2 

AT 
V 

V CT 

V(X),V(Y),V(Z) velocity components along X,Y,Z coordinate axes 

AV(X),AV(Y),AV(Z) velocity increment components 

AV^AVf initial velocity impulse; final velocity impulse 
total velocity impulse 

coordinates fixed to interceptor initial orbit (inertial) 
coordinates fixed to interceptor transfer orbit (inertial) 
x",y",z" coordinates fixed to target vehicle (rotating) 


AV 

x,y,z 


xW 
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X,Y,Z coordinates fixed to target orbit (inertial) 

a, 5,0 auxiliary quantities 

y constant determining hyperbolic or elliptic computation 

A0 transfer angle 

9 true anomaly of interceptor in transfer orbit 

69 increment of 

6t increment in transfer time 

6t w increment in wait time 

5v increment in true anomaly of interceptor in initial orbit 

I u gravitational constant, ft3/sec2 

v true anomaly of interceptor in initial orbit 

£ radial distance to interceptor in initial orbit 

£j(x),|^( y),ij(z) velocity components of interceptor immediately before initial impulse 

p radial distance to interceptor in target orbit 

p(x'),p(y r ),p(z’) velocity components of interceptor in transfer orbit 
(j> true anomaly of target 

\p auxiliary angle 

u> longitude of periapsis measured from node in plane in question 

ft longitude of ascending node measured in target plane 


Subscripts : 

c circular orbits 


f 


final time t = t 


f 


I 


initial time t = t ^ unless specified differently 
initial interceptor orbit 
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I 


k arbitrary index 

o observation time t = t Q 

p parabolic orbit 

T target orbit 

w wait time 

1,2 used to distinguish between the two terminal vectors in transfer orbit 

iterations 

1,2,3 used to indicate vector components along indicated Cartesian coordinates 

1,2, 3, 4 indicates velocity components immediately before and after initial (1,2) and 

final (3,4) impulses 

X,Y,Z used to indicate vector components along X,Y,Z coordinate axes 

Nonsubscripted Keplerian orbital parameters refer to the transfer orbit. Dots over 
symbols denote derivatives with respect to time; a caret (") over a symbol denotes 
vector quantities. 


DISCUSSION OF ANALYSIS 


The coordinate systems are chosen with the target orbit plane as the reference 
plane, the periapsis of the target defining the reference position vector in space. Fig- 
ure 1 is presented to emphasize the 
orbital trace of target parameters involved in describing 

the target orbit and the target vehi- 
cle. The target vehicle radius is 
denoted by r and its true anom- 
aly by (p . Input is given to define 
the geometry of the target orbit 
e,p, arp and the true anomaly of 
the target vehicle (p 0 at time 
t = t Q . The coordinate system 
X, Y, and Z is defined with X 
piercing the target orbit periapse, 

Z pointing along the positive rota- 
tion vector of the target, and Y 
completing a right-handed triad. 
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Unless otherwise noted, the quantities in this study are all nondimensional for gen- 
erality. This condition allows the deletion of two parameters p, and a^, from the per- 
tinent equations. The semi-major axis a T is not needed when dimensionless quantities 
are considered inasmuch as it is used for the normalizing and is only necessary when 
certain dimensional information is required. Neither is the central-body gravitational 
constant p required for dimensionless studies but must be used if the times are 
required to be dimensional (minutes, days). 


The parameters required to describe the initial interceptor plane are shown in fig- 
ure 2. These parameters are the usual Keplerian elements describing the initial plane 
and the periapse position of the interceptor orbit S2j, Wj, ij and the geometrical ele- 
ments of the interceptor ellipse ej, aj. The position vector of the interceptor in the 
initial interceptor orbit is given 
by the magnitude £ and the true 
anomaly v. The coordinates 
x, y, and z are also defined 
in this system; x being 
directed through the periapse, 
z perpendicular to the orbital 
plane (see fig. 2), and y com- 
pleting a right-handed triad. 

It is not necessary to 
specify the six parameters for 
the interceptor position in the 



Keplerian form (fij 


Ut 


ej, aj, and u Q \ as the program 


provides an alternate form for 
this input. The position and 
velocity (x£, y£, z”, x 


Figure 2.- Plane of initial interceptor orbit and associated quantities. 


Plane of interceptor 


y'A> 


and 




o’ J o’ 
of the interceptor rela- 


tive to the target at time t = t 0 
may be used in place of the 
Keplerian parameters. (See 
fig. 5.) A derivation of the rela- 
tionship between these param- 
eters is given in appendix A. 

The transfer orbit plane is 
specified entirely by the geom- 
etry of the initial and final posi- 
tion vectors of the interceptor. 
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(See fig. 3.) The time for the initial impulse is denoted by t = tj and the final impulse 
t = tj. The quantities <p Q and v Q are the true anomalies of the target and interceptor 
at some time t = t Q but the initial impulse may not come until some time later tj; 
thus, another quantity is introduced, the wait time At^ where = tj - t Q . It is then 
clear that the transfer time At is At = tj - tj. The initial position vector |j is deter- 


mined by the orbital quantities and the wait time ej, aj, fij, ij, cuj. 


v Q , and 


AW 


kZ 



Figure 4.- Quantities required for Keplerian input to program (DMAS = 1.0). 


or the relative coordinate data 


*•0’ 


y'o> 


z o’ 




yB> 


z”, and 


At w . To determine the final posi- 
tion vector ? p the position vec- 
tor of the target at t = tj, the 
quantities e T , <p Q , and At w + At 
are needed. The total time from 


to 


is the sum of the wait 


time and the transfer time. 

The initial position vector 
|j is determined by first com- 
puting the eccentric anomaly Ej q 
at t = t Q of the interceptor in 
its initial orbit from the equation 


tan 


E 


Io 


1 - e T V n 

rr^ tan T « 


AZ(K) 



Figure 5.- Quantities required for relative coordinate input to program 
(DMAS = 2.0). 


The mean anomaly at t = t Q is 
found directly from Kepler's 
equation as 

M Io = E Io - e I sin E Io ( 2 ) 

and at t = t j the mean anomaly 
will be Mjj where 

Mi i .M Io+ 2,rgAt w (3) 

The eccentric anomaly E^ at 
t = t j is found from iteration of 
the Kepler equation as described 
in reference 5. 
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(4) 


M n = E Ii + ej sin E n 

After obtaining from equation (4), the true anomaly may be found from 

E T 


v . n h, t . 

*“T 


After the true anomaly of the interceptor v j preceding the first impulse is found, 
the radius vector magnitude is obtained from the well-known solution 


5i = r+ 


a l(* " e I 2 ) 


ej cos 


(5) 


The rates may be easily computed from the relations 

i[ a i(i - e i 2 ) 
1 ' («i ) 2 

and 


( 6 ) 


k- 


ej sin ia 

W 1 - e i 2 ) 


(7) 


The initial vector 


written in the x,y,z coordinate system (fig. 2) becomes 


li = 5j(cos 


iai + sin 


-if) 


( 8 ) 


The final position vector is determined from the initial conditions e,p and <p Q 
along with the elapsed time At w + At in the same manner as the initial position vector. 
These computations are shown in appendix B. By repeating these steps for the target 
vehicle at t = tj the true anomaly may be obtained and may be used to find the 
final radius vector written in the X,Y,Z system as 


= r^cos <pf I + sin cp^j 


where 


r f _ 


e 2 
e T 


1 + erp cos 


(9) 


( 10 ) 
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Also 


0f — 


r f " 


^1 - e,p2 

(id 

w 2 


^ sin. (j!)| 

(12) 


1 - e r 


The interceptor velocity immediately preceding the initial impulse at t = t^ is 
obtained from 

ii = (lj cos ia - Sfi sin v^ji + ^ sin ia + ^ia cos (13) 

and after the final impulse the interceptor velocity must be the same as that of the target 
or 


?f = (r^ cos - r f 0£ sin <p^j I + ^r f sin + r f <£j cos cp^jJ (14) 

The transfer orbital plane and some of the properties of the transfer orbit are found 
by noting that the interceptor must leave vector at time t and arrive at vector r f 
at time tj. Since and r^ are known, the transfer plane properties and the transfer 
angle Ad may be determined. 

In equations (8) and (9), was specified in terms of its components in the x,y,z 
coordinate system whereas r f was specified in the X,Y,Z system. In order to manip- 
ulate with these vectors, their components must be referenced to the same coordinate 
system. The target coordinate system X,Y,Z is used here so that the vector must 
be transformed. This transformation can be accomplished by the common Euler angle 
matrix. The elements of this matrix are dependent upon the angles fij, c<jj, and ij and 
are given in appendix B. The matrix will here be simply noted as Ja^j . Then f ^ is 
written with components in X,Y,Z as 

= [ a ij]#i>k) (15) 


or letting 


li(i,J,K) = ^(X)I + £.(Y)J + ^(Z)K 


(16) 
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where ^(X), ^(Y), ^(Z) are functions of ^ , ia, u^, S^, and ij and are found 

from equation (15). These relations are written out in appendix B (eqs. (B21)). 

With both and r f in the X,Y,Z coordinate system, the transfer angle A 0 
and the transfer -orbit plane inclination i may be derived from the vector identities 

ii * ?f = £i r f cos Ad 
li X i-f = ^rj sin A0k' 


where 

A ^ A ^ 

k r = -sin i sin cp^L + sin i cos cp^J + cos iK 

A complete derivation of Ad and i is given in appendix B. 

The velocity components of equation (13) may also be transformed to the X,Y,Z 
system 

ii(U,K) = [a^KlJ.k) (17) 


which gives expressions for ^(X), ^(Y), and ^(Z) as shown in appendix B (eqs. (B44)). 

The problem now becomes one of finding an arc of a conic section which will pass 
through the initial vector and the final vector r ^ and which also has the property 
that a body traversing this arc will do so in the desired transfer time At. For future 
convenience, the following convention is defined: Let p^ be the minimum of r^ and 
and pg the maximum; then the transfer orbit must be the conic section passing 
through the radii p^ and p 2 , separated by an angle Ad, in the time At. 

If an angle 0j of periapse to pj is guessed, the eccentricity for the conic section 
which has radii p ^ and p 2 separated by Ad may be found from 


where 


p i = I7 


e cos 0-. 


P 2 "17 


e cos 0c 


(18a) 


$2 = + AQ 


(18b) 
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and 


p = a(l - e^) 


(18c) 


from which p is eliminated to find 


e 


p 2 ~ P 1 

Pj cos - p 2 cos 0 2 


The eccentric anomalies may also be computed from 


tan ^ = \ 


1 - 


1 + e 


e tan^ 


(19) 


and the mean anomalies from 


M = E - e sin E 

The time of transfer for this particular choice of 6 ^ is found from the mean anomalies 
as 

AT = £ AM 

2ir 

The semi- major axis a is determined from e and 0^ or 0g in equation (18c). 

The computed transfer time AT will not, in general, correspond to the desired 
transfer time At; therefore, it will be necessary to assume a new value for 0j and 
continue with this process until the desired agreement is obtained. 

It is possible to find certain regions from which to choose 0^ and thereby shorten 
the iteration considerably. It is easy to see in which conic section the transfer orbit must 
be by computing the time required to transfer by a parabolic orbit. This time may be 
determined without iteration by setting e = 1.0 and following the procedure outlined in 
appendix B. If this parabolic time AT p is less than the desired transfer time At, the 
orbit must be elliptic (e < 1.0) whereas if the parabolic time ATp is greater than the 
time At, the orbit must be hyperbolic (e > 1.0). The procedure is to find the type of 
orbit, hyperbolic or elliptic, and then set 0^ equal to the parabolic anomaly 0 p plus or 
minus some increment 60 j so that 0^ lies in the hyperbolic or elliptic region, which- 
ever is appropriate. The size of this increment 60 j is important in convergence of the 
iteration. It is governed by two other limiting values, one each for the elliptic and hyper- 
bolic regions. Then a straightforward incrementation process is used directly until the 
correct value of 0 ^ is found. A complete description of this process is found in 
appendix B. 
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After the properties e, a, 0^, A0, p, £2, u>, and i of the transfer orbit are 
found, the next step is to compute the velocity increments required for the rendezvous. 
The velocities before the initial impulse and after the final impulse are given in equa- 
tions (17) and (14). The velocities after the initial impulse and before the final impulse 

A. A X 

are described in the transfer orbit. If i T , j', and k T are the unit vectors along the 
x* ,y',z' coordinate axes previously defined in the transfer orbit system, these velocities 
are written as 

P* = cos e i - P i e i sin 0^i' + ^ sin 9 i + P i 9 i cos 9^]' 

Pj = (pf cos 0j - p 0j sin 0jji' + ^p f sin 0j + p f 0j cos 0 f jj' 

These velocities are then transformed to the X,Y,Z system for ease of manipulation. 
This transformation is done by a matrix identical functionally to [ay] but with 

the transfer orbit angular parameters £2, a>, and i replacing the parameters for the 
interceptor initial orbit £2j, Wj, and ij. After these transformations the velocity incre- 
ments are easily found from 


Av i - h - k 


AV f = r f - p f 

and the total velocity increment required becomes 


AV = 




2 


A thorough discussion of this method and the pertinent mathematics of the problem 
is found in appendix B. 


COMPUTER INPUT AND OUTPUT 

The problem was programed in the FORTRAN IV language for the IBM 7094 com- 
puter installation at Langley Research Center. Certain options were incorporated in the 
input of the program and are described. A more thorough description of the entire pro- 
gram is found in appendix C. It was felt that the present writeup is sufficient for most 
uses and if changes are desired, they are easily incorporated through the FORTRAN IV 
language. 
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The program automatically increments the anomaly v Q of the interceptor at 
t = t 0 , the wait time in orbit At w , and the transfer time At. The incrementation begins 
with At or At^, depending on the value of the option control variable GUIDE in the 
input. The second parameter incremented is At w or At, whichever is not incremented 
first, and the third incrementation is v Q . 

An option is also available for using the input and output times At and At w as 
dimensional quantities. The input quantity OPTION determines whether the quantities 
At and Atw are dimensionless or dimensioned as minutes or as days. Because of the 
nondimensional parameters used in the program, the input ju. and a,p may be any con- 
venient values as long as the time is dimensionless and Prp and are not desired; 

otherwise, \i and a^, must be properly dimensioned. 

The input quantity DATAS determines whether the initial interceptor orbit data are 
Keplerian or relative Cartesian. 

Other necessary input quantities are h, the increment constant for the itera- 
tion (in the examples solved a value of 3 was used), 6t and 6t w , the transfer and wait 
time incrementation for the time grid. Also, input is the quantity r min which is the 
minimum radius the interceptor may take as it travels the transfer arc. This value does 
not affect the computations in any way; however, if the radial distance of the transfer arc 
falls below the value r m ^ n , an asterisk is printed out at the right-hand side of the output 
sheet. 

The variables output in the program described herein were those which were con- 
sidered to be of general use. The output is printed as shown in tables I, n, and in which 
are output data for the three examples following this discussion. The first block of data 
is a reprint of the input where the following correlation between symbols may be observed 


Symbol 

Corresponding value 

ET 

e T 

El 

e I 

AI 

a I 

OMEGA I 


OMI 

"I 

n 

h 

PHIO 

*0 

NUOl 

v Q (first) 

NUOL 

v Q (last) 

MU 


AT 

^ (ft) 

TPER 

P T (sec) 

VCT 

V CT (ft/sec) 

DTTI 

At (first) 

DTTF 

At (last) 

DTWI 

Aty, (first) 

DTWF 

AT W (last) 

RMIN 

r min 

H 

h 


In the case that the input is in the form of relative 
coordinate data, these Keplerian quantities for the initial 
interceptor orbit are computed from the relative coordinate 
data as has been explained. An example of this printout is 
shown in table H where a fourth row of variables is noted and 
XO, YO, and ZO correspond to Xq, y^, and Zq and 
XODT, YODT, and ZODT correspond to Xq> Yo> anc * zJJ. 

After this output which records the data for the terminal 
orbits, the data are computed for the rendezvous problem and 
the following description of the printout nomenclature should 
be noted: 
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Symbol 

Corresponding value 

NUO 

u o 

DELTTT 

At 

DELTTW 

at w 

PHIF 

(j)f 

NUI 

v i 

E 

e 

A 

a 

I 

i 

THI 

«i 

DELTTH 

Ae 

RHOM 

p min 

DV1X 

Component of AV^ along I 

DV1Y 

Component of AV^ along J 

DV1Z 

Component of AV^ along K 

DV2X 

Component of AVj along I 

DV2Y 

Component of AV| along J 

DV2Z 

Component of AV^ along K 

DELV1 

Magnitude of AV^ 

DELV2 

Magnitude of AVj 

DELV 

Magnitude of AV 


EXAMPLES 

Example 1 — Interplanetary Transfer Between Massless Planets 

To compute velocity increments required for interplanetary trajectories the input 
for the planet ephemeris is needed. This input can be found from reference 6. As an 
example, consider a trip from Earth to Mars which is to be accomplished by a two-impulse 
orbit. As the observation time t Q , January 4, 1964 (Julian date, 243-8398.5) is selected. 
From the 1964 ephemeris tables (pp. 18, 50, 172, 176), data for this problem are found 
which are referenced to the ecliptic plane. From the geometry of the ecliptic coordinate 
system (see refs. 1 and 2) and the definitions presented here, the following program input 
is determined: 

e T = 0.093372 v Q = 0.37° 

a T = 7.4737 X 10 11 Oj = 253.88° 
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0 o = 324.4° cDj = 233.02° 

ej = 0.0167242 ij = 1.850° 

aj = 0.656301 

In order to be notified in the output whether the transfer arc is less than 0.6 astro- 
nomical units, let r . = 0.6a T . Also let 

’ mm I 

p = 4.679 x 10 21 ft3/sec2 


h = 3.0 


then. 


r - 
mm 


0.393781 


Since v Q is not to be incremented (because it is fixed by the physics of the solar 
system), let v Q (first), v Q (last) be v Q as shown and the increment 6v Q = 0. Let it be 
desirable to collect data in increments of 20 days with transfer times from 160 to 260 days 
and waiting times of 0 to 40 days, and then use 


At(first) = 160.0 At(last) = 260.0 6t = 20.0 
At w (first) = 0.0 At w (last) =40.0 6t w = 20.0 

If the data are to be analyzed with At as the primary independent variable and At^ as 
the parameter, At is incremented first. Therefore, set GUIDE at 2.0. The times are 
to be input and output in days so that OPTION is 3.0, and since Keplerian input is used, 
set DAT AS as 1.0. Some output data are given in table I with time in days and all other 
quantities except p, a^,, Prp, and V CT are nondimensional. 

Example 2 - Earth- Orbit Rendezvous 

There are many approximate schemes which have been derived and investigated for 
determining the velocity increments required to rendezvous between similar orbits. Such 
investigations have been mainly concerned with circular orbits or first-order representa- 
tions of elliptic orbits (near circular orbits). The simplest case of rendezvous with a 
circular target orbit (as done by Clohessy and Wiltshire, ref. 7) gives a closed-form solu- 
tion to the linearized equations of motion for the velocity increment required to rendezvous 
in a given amount of time. This linearized solution is accurate for the terminal phases of 
a general rendezvous problem but the accuracy is soon lost for large initial interceptor- 
target ranges. It should be noticed that even in the docking phase, this solution may not be 
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too accurate if the target orbit is highly eccentric. Further attempts to increase the 
applicability of the Clohessy-Wiltshire results usually require an iteration scheme 
(Anthony and Sasaki, ref. 8) or some other process to get a solution of the velocity incre- 
ment required to rendezvous. 

It is noted that the solution of the exact equations presented herein is not too 
involved once a computer program is available, and, of course, offers the advantage of 
giving correct results for both large eccentricities and inclinations. 

As an example of the use of this program for earth-orbit rendezvous, assume that 
the following elements are used for the target orbit: 


Also let 


e T = 0.0234 
a T = 2.248 X 10 7 ft 


0 O = 0.0° 




r . 
min 


0.954 


1.408 x 10 16 ft3/sec2 


h = 3.0 


The interceptor position and velocity are given in relative coordinates : 


x£ = -0.01692 

y” = 0.0376 
J o 

z" = 0.0 


x” = -0.00376 

y" = 0.1526 
J o 

z” = 0.0 


Note that ep a^, v Q , ip £2p and o)j are computed from these relations. Let the time 
be input and output in minutes: 

At(first) = 2.0 At(last) = 40.0 6t = 10.0 

At w (first) = 0.0 At w (last) = 30.0 6t w = 10.0 

If a plot of At^ as abscissa and At as a constant parameter is desired, set 

GUIDE at 1.0. Since time is in minutes, OPTION is 2.0, and since relative coor- 

Output 

dinate data are used, set DATAS at 2.0. 

The output of the computer program for this input is given in table n. 
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Example 3 - Transfer Orbits 

The problem of orbital transfer between two arbitrary elliptic orbits may also be 
accomplished by this computational procedure. The use of this method in orbital transfer 
problems is indicated in this example where the minimum two-impulse transfer orbit is 
desired between two fixed terminal orbits. The transfer orbit requiring minimum velocity 
increment may be found by proper interpretation of the data and the use of a plotter. The 
program as presented here does not have a gradient or other optimization scheme to find 
the minimum velocity transfer directly. However, an example is shown of one way of 
obtaining the optimum (two- impulse) orbital transfer. 

For the terminal orbits the following input data are used: 

e^, = 0.5 arp = 1.0 (arbitrary) 

= 90° e r = 0.2 

ij = 30° a r = 0.9 

Wj = -90° 

Also let 

r . = 0.45 
mm 

h = 3.0 

iU = 1.0 (arbitrary) 

The true anomaly v Q of the interceptor at t = t Q is chosen in the following range: 

^ 0 (first) = 0.0 
v 0 (last) = 350° 

5v Q = 10° 

which is incremented automatically by the program. To obtain all the possible transfer 
arcs, it is necessary to increment either or v Q , but not both. As v Q is chosen 

to be the incremented parameter, the wait time is set equal to zero. 

The true anomaly 4> Q is chosen, and when the transfer time At is varied, a set 
of curves for the variation of AV with At with v as the parameter are obtained for 
a specified value of (p Q . The program does not increment <p Q automatically as it does 
for At, At^, and v Q . However, it always returns to read the first data card of a new 
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set after running through one complete set of data. If there are no more sets of data, the 
computation ends; if there is another complete data set, it starts anew. To increment 
c p , it is necessary to reproduce all data cards for every new value of <fi Q and to place 
these one behind the other. (See appendix C.) 


As set out above, each complete data set gives one complete set of curves for a 
given <p Q with v Q as the parameter. (See sketch 1.) 

As an aid, sketch 1 has been drawn 
with <p£ as the abscissa value rather 
than At and is possible because of the 
functional relation between At and < 
or 


<Pf = <t>o + F ( At ) 

since At^, = 0. 

The transfer time At is read in 
as follows : 



At(first) =0.09 




At(last) = 1.008 
6t = 0.009 


Sketch 1 


These quantities are dimensionless ratios to the target orbit period. The waiting time is 
input as 

At w (first) = At w (last) = 6t w = 0 


Some sample computer output is shown in table in for <j> 0 = 0.0°. 

From a Beckman automatic plotter, the characteristics of one of the lower AV 
transfer orbits are found to be 


v i = 55° 
c t> t = 193.58° 
At = 0.64 
e = 0.4361 
a = 1.086 
6 i = 71.96° 


A0 = 132.47° 
i = -22.88° 


Pmin = °- 78 


AVj = 0.0808 

AV f = 0.2563 
AV = 0.3371 
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CONCLUDING REMARKS 


This report presents a technique for solving three-dimensional orbital problems in 
a straightforward manner using the exact solutions or Kepler solutions to the equations 
of motion. Basically the method developed is an iteration on the Kepler equation using the 
true anomaly as the iteration parameter and the mean anomaly or transfer time compared 
with a prespecified transfer time as the stopping criteria. To aid in the choice of the true 
anomaly to begin the iteration, certain boundaries are devised within which the solution 
must lie. The iteration is performed directly on the Kepler equations and no derivatives 
are necessary. This method works very well and the computation time compares favorably 
with other methods, the typical run times being about 0.003 minute per transfer. The 
generality of the program format presented allows rapid computations, with simple engi- 
neering input parameters, of interplanetary rendezvous or near -planet rendezvous cases. 
For instance the program has proven useful in studies concerning fuel requirements for 
abort missions during lunar letdown of the lunar excursion module and command module. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., February 14, 1966. 
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OBTAINING KEPLERIAN ORBITAL ELEMENTS FROM 
RELATIVE COORDINATE DATA 

In many problems concerning rendezvous, the initial conditions are given in terms 
of the relative position and velocity of the interceptor with respect to the target. If the 
Keplerian equations are used for the exact solution of the motion, the Keplerian elements 
in terms of the relative coordinates must be obtained. Let be the relative posi- 

tion vector (see fig. 5) of the interceptor with respect to the x",y",z" coordinate system 
attached to the target; x" points away from the gravitational center, y" is in the posi- 
tive direction of angular motion, and z" is along Z normal to the x"y" plane. 

Suppose that the following are given in dimensionless form: e T , <p Q , x£, y”, z’^, 
x" y'^, and z'^. It is desired to find the Keplerian constants ej, aj, v Q , S2j, coj, 
and ij. By using the formal application of vector analysis, let i,j,k be unit vectors 
along x", y", and z", and obtain the radius vector 

£ = (r + x")i + y"j + z"k (Al) 

for the interceptor, and since the rotational rate of the (x”,y",z") coordinate system is 
(pk, the velocity is 

| = (r + x')i + y'j + z'k + <£k x f (A2) 

where 

k X £ = (r + x")j - y"i (A3) 


The velocity of the interceptor becomes : 

V = | = (r + x" - + (y” + <^>(r + x"))j + z"k 

The fundamental relations of orbital mechanics show that 


P T = 1 - e T ^ 

Pip 


r = 


1 + e>p cos 4> 


. e T sin <p 


PT 


<P 




(A4) 


(A 5) 
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From these expressions the velocity of the interceptor in its initial orbit can be written 
in terms of the desired quantities along the axes of the moving x ,T ,y",z" coordinate 
system. Quantities considered here are at the time t = 0 and are subscripted accord- 
ingly. The velocity components at t = t Q along the x", y", and z" axes will be 
denoted by V Q j, V^, and V Q g, respectively; thus, the interceptor velocity becomes 

v = V + V o2i + V o3* < A6 > 


By using equations (A4), (A5), and (A6), the components are written in terms of known 
quantities as 


~\ 



The energy equation for the interceptor at time t Q may be written (nondimensionally) as 

V 0 2 -2-=-± (A8) 

^o a I 

where £ Q , V Q are magnitudes of the defined vectors £, V at time t Q or 


«o - + x o) 2 + K) 2 + Kf 

v o » fol 2 + ^7+^? 

By using equation (A8), solve for the semi-major axis aj 



(A9) 

(A10) 


(All) 


To determine the eccentricity of the initial interceptor orbit, write the angular 
momentum at time t Q of the interceptor, 


H o^o* V o 


(A 12) 
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where | and V Q are dimensionless quantities. If 

«o = H ol‘ + H o2? + H o 3 £ 
performing the indicated operations yields 

H ol - y’o V o3 - z o V o2 

H o2 = z o V ol - ( r o + ?o) V o3 ) 

H o3 “ ( r o + x o) V o2 “ y’o V ol 

and 

H o ■ f H ol ) 2 + ( H o 2) 2 + ( H os) 2 


(A13) 


(A14) 


In the nondimensional form, the angular momentum and semi-latus rectum may easily be 
shown to be related as 



(A15) 


and therefore the well-known relationship pj = ajfl - ej2j gives with equation (A15): 


K, k 



e I = 




(A16) 


Sketch 2 and figure 5 are given as aids in 
describing how the angular measures ip Oj, 
and u>j are determined. The unit vector N 
directed along the line of nodes of the initial 
interceptor orbit and pointing toward the 
^ ascending node is defined by the relation 

kXH Q = |h q sin ij|N (A17) 

The inclination ij is also described by the 
relation 


k • H Q = H Q cos ij (A 18) 
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and equation (A18) gives 


i T = cos 


_lfo3 
H„ 


From figure 5, it is geometrically evident that 

I = i cos (j) Q - j sin 4> Q 
and also from the definition of the angle Oj 


cos — I * N 




k sin S2j = I X N 


and expanding the vector manipulation of equation (A17) gives 


n = . j - _ H ° 2 , i 


From the well-known relationship 


obtain 


H Q sin ij H q sin ij 


{ h 

° 1 + ej cos v Q 


1 y I i 
COS - 1 

° e I^o 


A basic orbital relationship which is of value here is 


Sin V n = IP ir 


o ej ^o 


and, as is the component of V Q along | Q , 


f = V • — 
? o ot 
? o 


(0 < i 


1 I < 


In order to find the angle u>p note that another geometric relation is 

j o -N = 4 0 cos (^+^ 0 ) 


7 r) (A19) 

(A20) 

(A21) 

(A22) 

(A 23) 

(A24) 

(A25) 

(A26) 

(A27) 
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and also 


H 


NX S 0 = ^ 0 sin ( w I + ^o)h 


Note also that v may be found from equations (A24) and (A25) as 


v = tan" J 


( Os "o < *4 

From equations (A20), (A21), and (A23), 

H o2 cos ^o + H ol sin K 


cos S2j = - 


H o sin h 


and the proper vector manipulations with equations (A20), (A22), and (A23) yield 

H ol cos 4>o “ H o 2 sin 4>o 


sin = 


H o sin h 


(A 28) 


(A29) 


(A30) 


(A31) 


From equations (A30) and (A31), 


= tan 


-1 


^sin SI 


I 


\cos S2t 


(0 ^ S2j < 2tt) (A3 2) 


Similarly, performing the indicated vector manipulations on equations (A27) and (A28) 
leads to the scalar equations 


cos la, + ,j = ^2 H 0l - S ol H o2 
1 1 ° j «o H o sin 4 


sin = 


So3 


?o sin h 

If the inclination ij is zero, S2j is arbitrary; thus, let 


S2j = 0 


Also, equations (A27) and (A28) yield 


cos 


h + v o) 


*ol cos ^o - Sq 2 sin *o 


(A3 3) 
(A34) 


(A3 5) 


(A36) 
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sin ( W I + v o) = 


gpl sin K + ^o2 cos K 


Whether ij is zero or not, equations (A33), (A34), or (A36), (A37) give 

~sin (qjj + iO 


Wj + = tan 


-1 


cos (coj + V Q ^j 


and equation (A38) with equation (A29) give the desired expression 


“I = ( W I + v o) ~ v o 


(A3 7) 


(A38) 


(A39) 


Equations (All), (A16), (A19), (A29), (A32), and (A39) give the information necessary 
to compute the desired Keplerian elements. 


26 


III II 


APPENDIX B 


MATHEMATICAL DESCRIPTION OF PROBLEM 


The following description is a logical flow of the problem as it is programed. 


Description of Initial and Final Properties 

If the initial and terminal states of an orbit referred to two other elliptic orbits, 
namely, the interceptor initial orbit and the target orbit, and also the time required to 
transfer from orbit to orbit are known, the velocity increments required to establish the 
transfer orbit and the elements of this transfer orbit may be computed. The initial and 
terminal states may be found from the input as follows. Suppose that the following data 
are given: ju, a^, e-p, <£ 0 , ej, aj, coj, Oj, ij, v Q , At w , and At. Recall that all 

quantities except p and a,p are dimensionless. 

Compute the semi-latus rectum of the initial interceptor orbit and the target orbit 
from 

P I = a l( X ‘ e I 2 ) (Bl) 


P T = 1 - e T 2 


(B2) 


The target orbital period in seconds is 



(B3) 


The eccentric anomaly of the target at t = t Q is found from e^, and <p Q by the 
defined relationship 



1 - e 
1 + e 



and then functionally as 


E To “ F ( e T’^o) 

For computational purposes, define the function F(a;,/3) as 


F(a,0) = 2m + 2 tan 


-lM i - oT 
\ 1 + a 


y 


tan ^ 
2 


I s *- and 4*-* 1 quadrants 


(B4) 


(B5a) 
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where 



and 

n = Integer ^ + v \ 

\ /lowest integer value 


(B5b) 


(B5c) 


y = 1.0 (a < 1.0)1 

> (B5d) 

y = -1.0 (a > 1.0 )J 

Note that allowance is made here for the possibility of hyperbolic orbits because F(o?,/3) 
will also be used for transfer orbits. 

The mean anomaly of the target orbit at t = t Q is found directly from Kepler's 
equation 


M>p 0 — ^tq (B6) 

The mean anomaly of the target at the final time t = t^ is 

Mrp£ = Mrp 0 + 27r^At + At^j (B7) 

The eccentric anomaly of the target at t = t^ is required so that the true anomaly 
at t = t f may ultimately be obtained. It is possible to solve for the eccentric anomaly by 
use of an expansion in terms of the infinite Bessels series (see ref. 9) but it was found to 
be more rapid to iterate the Kepler equation by using the method of differential correction 
as shown in reference 5 or 9. The steps in the iteration procedure are as follows: 

(1) Estimate the initial value E^ from truncated series solution of the Kepler 
equation 

E 1 = M o + e s * n M o + 2 Sin ^ M o ( M ° = 6 = 6T ) (B8) 


(2) Compute the following sequence: 
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i = 1 
M, = E, - e, 


T sin Ej 


l 1 

AM^ = Mrj,j - 

AEj = AMj j(l - e T cos E^ 


E . j —— E . -|- ^ ij . 
1+1 1 1 


Compare 


E i+1 : 


i = i + 1 


|E i+ l * E. 
\E i+1 = Ei 


which gives after the completed iteration: 


E 


Tf 


= E. . 
l+l 


(B9) 


At t = tj the eccentric anomaly of the interceptor E^ is obtained from v Q , ej, 
and At w in the same manner as E Tf . The expressions are rewritten here for com- 
pleteness with eccentric anomaly of interceptor at t = t Q as 



where F^ej,^ is the function previously defined (eq. (B5)) with a = ej and /3 = v Q . 
Also 


M Io - E Io - e I sln ( E Io) 

At t = t^ the interceptor mean anomaly is found from 

T 


M Ii = M Io + 2n ' 


3 At w 
aT" 3 


The iteration steps are shown in equations (B8) and (B9) with the following 
replacements: e T by ep Mrpj by Mj.; Erpf by Ejp 

With this information the true anomalies at the initial and final times of the intercep- 
tor in its initial orbit and the target in its orbit may be computed. The routine given by 
equations (B4) and (B5) will give these values if -a is replaced by the corresponding e. 
Thus, the true anomaly of the interceptor in its initial orbit at t = tj is, a and (3 in 
F being replaced by the parameters -ej and Ejj 

v i = F(- ei ,E n ) (BIO) 


and the target in its plane at t = tj is 
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cpf - F(-e T ,E Tf ) 


(Bll) 


With these quantities it is possible to compute the properties ja , r^, r^, 

and cpf which will be required for computation of the velocity increments. The initial 
radius vector is 


^ 1 + 


Pl 


ej cos Pj 


The rate of change of true anomaly at t = t^ is 


^ = ^ (si ) 2 


Also the rate of change of the radial distance is 

*i" 


e I sin *i 


fl 

For the final time t = t^, the properties of the target orbit are 

Pt 


r, = 


f "i 


+ e T cos 


(B12) 


(B13) 


(B14) 


(B15) 


- 



(B16) 


erp sin <p£ 

^T 


(B17) 


If an interceptor is to transfer (or rendezvous) from one orbit to another, the two 
orbits being determined by this information, it must begin at t = t^ in the initial orbit, 
change velocity (instantaneously) to the transfer orbit at t = tj, and travel until t = tj in 
the transfer orbit. Then it instantaneously changes its velocity to satisfy the terminal or 
target orbit properties. Enough information is now available to compute the initial and 
final terminal velocities but further discussion of this procedure will be delayed until the 
initial and final transfer orbit properties have been computed. 


Computation of Transfer Orbit Properties 

A* A 

Information giving initial and final terminal vectors in space r^ is available, 
and it is now necessary to determine the transfer conic which connects these two vectors 
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in the specified transfer time At. Reference 3 shows that the solution to this problem 
is unique. However, because of the implicit nature of the Kepler problem, it is necessary 
to perform an iteration in order to obtain the solution. The method used here is a 
straightforward iteration of the Keplerian orbital equations as was also described in 
reference 3. 

The plane of the transfer orbit is specified by the plane containing the initial and 
final position vectors and the center of attraction. The transfer plane inclination i and 

A 

the transfer angle A0 are determined with the aid of vector representation. Let I, 

J, and K be unit vectors along X, Y, and Z axes; let i, J, and k be unit vectors 
along x, y, and z axes as shown in figures 1, 2, and 3. Let Ja^TJ be the x to X 
transformation matrix so that a general vector A transforms as: 

A(IJK) = ja i;j JA(ijk) 

The initial vector may be written as 

h = ^i( cos V + sin ( B18 ) 

By letting com P onen t s of fj along the X,Y,Z axes, the following 

relations are obtained: 

k = *il* + *12* + ^i3* 

fj = r^cos <p £ I + sin 

In order to examine the vectors in the I,J,K system, a transformation of components is 
required; thus, 

^(I.J.K) = ^^(ijk) (B19) 

where 



COS 

cos - cos ij sin sin uij 

- sin ct>j cos - cos ij sin cos coj 

sin Oj sin ij 

N- 

cos 

sin + cos i^ cos sin 

- sin cdj sin + cos ij cos cos 

- cos sin ij 


sin ij sin oij 

sin ij cos 

cos i T 

1 J 


From equations (B18), (B19), and (B20), 
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^il = ^ip os cos (^i + - cos ij sin £2j sin^ + a>jjj 

^i2 = ^ij® in % cos^j + + cos ij cos sin^i/j + cojjj 

*13 = Ijjsin ij sin^ + ^j] 

To determine A 0 and i, note the vector identities: 

li • r f = ^r f cos A0 
^i x r f = ^r f sin A 0k' 

Note that the unit vector k' along the z' -axis is written as 

k' = I sin i sin 0 - J sin i cos £2 + K cos i 


and since 

Q + ij = (pf 

k' = - sin i sin I + sin i cos <p^J + cos iK 
From equations (B18), (B19), (B23), and (B24) find 


cos A0 = 


lil cos 0 f + £ i2 sin <p f 


k 


sin A 0 sin i = 


^i3 

*i 


tan i = 


*i3 


li! sin cj) f - £ i2 cos 0 f 


Combining these equations and solving for i yields the following expression: 

*i3 


i = arc tan 


lil sin 0 f - ^ i2 cos 0 f 




(B21) 


(B22) 

(B23) 

(B24) 


(B25) 
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and also 

cos Ae . Ill cos 0 f + i 12 sm 0 f 

«1 


S n sin <p, - 4 l2 cos <p, 
sin Ad = — 1 4 1 

COS 1 


(B26) 

(B27) 


Ad = arc tan 


sin Ad 
cos Ad 


The nodal angle 


52 may be expressed by 


52 = <pf + ir 
57 = <p£ 


(0 s a M 2tt) (B28) 


(i > 0) 
(i < 0) 


(B29) 


An inherent symmetry in the equations exists so that it makes no difference in the itera- 
tion whether the transfer from ^ to r f or r f to is considered, and advantage of 
this symmetry is taken by considering transfers only from the shortest of ff, ?f to the 
longest. Therefore, let pj be the minimum of rf and p 2 be the maximum of 

r f . (K r f = |f, Pf = P2 = P = r f = £i-) 

The chord c is given by 


c 


= \Jpi 2 + p 2 2 


2 PjPg cos Ad 


A geometric picture is shown in the following sketch: 



Sketch 3.- The quantities Pp P 2 , Ad, and c. 


(B30) 


By assuming a periapse angle of 6^ in the transfer plane, an orbit between p^ 
and pg can be determined as shown in the following sketch: 
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Sketch 4.- Transfer orbit for corresponding 0^. 


For each properly chosen angle 0p there is a conic section which passes through 
PpP 2 and an associated time of passage TpTg. Hence, for a given 0j, a given trans- 
fer time Tg - T^ is generated. In general, this time will not be the desired transfer 
time At, and thus an iteration is necessary. The necessary calculations for the deter- 
mination of the orbit by this method are now developed. Define 


e 2 = e 1 + Ad 

By eliminating p from the well-known expressions 

= P ^ 

1 + e cos 6 i 

p » 

1 + e cos 

J 

the unknown eccentricity of the transfer orbit is found from 


(B31) 


(B32) 


p 2 ~ P 1 

Pj COS - P 2 COS 0 2 


(B33) 


and the semilatus rectum 


p = p-j^l + e cos 0^ 
The semi-major axis is computed as 


(B34) 


a = — (B35) 
1 - e 2 

All the orbital parameters may be fixed by choice of the periapse angle 0j. 

A distinctive value of 0^ is the periapse angle 0 p for parabolic transfer from p^ 
to p 2 and the associated parabolic transfer time AT p = T 2 - Tj. This value of 0^ is 
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easily found since e = 1.0 for parabolic orbits. Letting 0^ = 0 p and using equa 
tion (B32) yields 


p^l + cos = p 2 (l + cos 0 2 ) = P 2 fl + cos^A0 + 


or 


(Pl - p 2 cos A0)cos 0 p + ( p 2 sin A0)sin 0 p = p 2 - Pj 


Then 0 is found by trigonometric identity 


sin^0p + i p'j = 


p 2 " Pi 


> 0 


where the angle \p is defined by 


sin = 


Pj - p 2 cos A0 


cos \p = 


p 2 sin A0 


* = tan-'fsJaLSii 
\cos l// 


2L < 
2 _ 


Let 


. -l/ p 2 " p l' 
01 = 3111 — 


and then note both values of 0 p in an interval of 2-n 


0 p = e p 


= a - \p 


minimum 


*P = 0 P 


— IT - (fit + \p) 


[maximum 


A geometrical description of if; is shown in the following sketch: 


(B36) 


(B37a) 



(B37b) 

(B37c) 
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Sketch 5.- Geometric interpretation of \p. 

The semi-latus rectum corresponding to 0p is 

Pp = Pl (i + cos e p ) 

The parabolic transfer time is then found as follows: The time to travel from periapse 
to p^ is found in reference 10 and is 

T i = ^-\jp 3 (i tan 3 + i tan 
Pi 2ir r P \6 2 2 2/ 


and from periapse to p 2 is 


i l — x/i „ + A0 

, = -r-\ p 3 4 tan 3 -2— 


0 + A P)\ 


so the parabolic transfer time is 


AT p - T p2 - T pl < B38 > 

This value may be used to compare the desired time At of transfer and to deter- 
mine whether the orbit is elliptic ^At > AT p j or hyperbolic (At < AT p j . 

There are certain definite regions in which 0^ must remain for a solution to exist. 
These regions are found by use of equation (B33) and the positive sign of e and p 2 - Pj. 
These equations imply 


Pj cos 0j > p 2 cos 0 2 
which leads to sin (o j + \j}j > 0 or 

i t - if/ > 6^ > - rjy 
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where 4 / is defined in equation (B37). It is seen that the lower limit for 9 ^ is - 1 // 
for all cases. The angle 9 ^ is greater than - 4 / as it is given by equation (B37). For 
an arbitrary e, 9 is replaced by 0, and equations (B37a) are written 

r 


0^ = sin" 


_i/ p 2 " P 1 


ec 


First quadrant 


- P 


From equations (B39) and (B37), the following inequalities are obtained: 


(B39) 


” e pU 


mimum 


e i < 0 p . . 

* minimum 


(e < 1) 
(e > 1) 


(B40) 


It may also be shown that the true maximum value for 9 ^ is less than n - 4 / and 
is the second parabolic solution 


0 = 0 
^ maximum p 


(B41) 


maximum 


Hence, the following regions for choice of the iterative periapse angles 9 ^ are 
determined from these relations: 

(1) If At < AT , then e > 1 (hyperbolic transfer) and -\J/ < 0. < 0 

r Ir 

(2) If At > AT , then e < 1 (elliptic transfer) and 9 < 9-, < ir - 2(\js + 0 \ 

Ir r A \ ir I 

(3) If At = AT p , then e = 1 and the solution is 0^ = 0 p 

A geometric interpretation may be given to these regions and an example is shown in 
the following sketch: 



Sketch 6.- Regions of possible choice for the periapse 9y 
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Elliptic transfer orbits .- In the case that At > ATp, the transfer orbit for rendezvous 
must be elliptic and hence 0 p < 0^ < 77 - 2ij/ - 0 p . The iteration scheme used here is a 
straightforward computation of the values Tg - T^ for values of 0^ starting with 0 p 
and continuing until Tg - Tj becomes equal to At. That is, the interval of 0^ is 
divided into h smaller intervals 60 where 


The first estimation of 


60 = 


0j is 


■n - 2(xf/ + 0 p j 



- 


p 

max P 

min 


h h 

0 + 60 and with this value compute 

r 

0 2 = + A0 


e 


p 2 ~ P 1 

p 2 cos e i - p 2 cos 0 2 



+ e cos 0 


1 ) 


a = 


P 

l- e 2 


and the eccentric anomalies 
to give: 


Ej and E 2 from the 
Ei = Fje^) 

E 2 = F ( e ^2) 


F(q!,/ 3) function defined previously 


and the elliptic time 


T 2 -Ti=ip 


E 0 - E 


j - e^sin Eg - sin Ejj 


This value of Tg - Tj is compared with At and if still too small, another increment 
60 is added to 0^ and the process repeated. When Tg - Tj is larger than At, then 
the interval 60 is halved and subtracted from the last value of 0j. This method of 
iteration is free from the singularities involved in methods using derivatives and is not 
significantly longer. The value of h used affects the time of iteration; however, an 
arbitrary value of h = 3 seems to produce sufficiently fast convergence. 
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Hyperbolic transfer orbits .- A value of At < AT p gives hyperbolic transfer orbits 
and hence -if/ < 6^ < 0 p . The method of iteration is identical to that for the elliptic case 
although the equations differ slightly. These equations appear as 

'P + 0 D d - (-ip) 

6 e = E = -2 

h h 

= 9 p - 59 

9 2 = 0^ + A 9 

The iteration begins at 0 p as in the elliptic case but now proceeds into the region of 
smaller 0^ 


e 


p 2 ~ p l 

Pj cos - p 2 cos 9 ^ 


p = p^l + e cos 0]J 


a 


P 

1 - e 


2 


(a < 0) 


The hyperbolic time func tions m ust now be used and equation (B5) is applicable if y is 
set equal to -1.0; thus, is obtained. Proceed with the usual hyperbolic "eccentric 

anomalies" which are the F(a,/3) functions with the preceding replacements. 


F i " F ’( e > 9 i) 

F 2 * F '( e ’ s 2) 

and the times from 0=0 to 0 = 0^ and 0 = 02 as 


T 


1 “ 


_ 1 _ 

2tt 



T 2 = ^f^^ tan F 2" 


L 
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Compare T 2 - Tj with At. If it is still larger, compute again with a new 
decreased from the last by 60. If T 2 - Tj is greater, increase 0j by and pro- 
ceed in this manner until the desired agreement is obtained. 

Special ca ses.- In the case that p^ = p 2 in either the elliptic or hyperbolic case, 
special computations are necessary. For instance, the iteration can no longer be accom- 
plished by incrementing 0^ as 0^ becomes a fixed value -\p. The problem is solved 
by iterating the eccentricity in a straightforward manner and is easily followed in the flow 
chart in appendix D. 

In the case that the time AT p happens to be equal to At, then there is no further 
iteration as the solution is the parabolic case. 

Computation of the properties.- After obtaining 0j, information on the rates of change 
is obtained. Also the quantities Pj and p 2 must be reassociated with r^ and 
The properties of the transfer orbit occurring at t = tp t = tj are defined as p- v Pp 
dp and 0j where 


Pi = Pi 

\k < r f 

Pi = p 2 

(«i =• r f 

Pf = P 2 

( r f =* k 

p f = P 1 

( r f •= k 

0i = 01 

[k < r f 

®l = " 0 2 

{k > r t 

0f = e 2 

( r f =• k 

II 

1 

H- 1 

( r f < k 


These relations follow directly from the geometry of the transfer orbit and the defini- 
tion that pj is the minimum of r^ and ^ and p 2 is the maximum of r^ and £p 
The rates of change at t = t p - t = t ^ in the plane of the transfer orbit may now be com- 
puted. Immediately after the impulse at t = tj, 





and immediately before the impulse at t = t^, 
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Pi = \fe esin0 f 


The transfer orbit plane is defined in exactly the same way as the initial interceptor 
orbit with the elements co, i, and 8\. The anomaly 8\ and inclination i were 
found previously. Once and a> have been determined, the transformation matrix 
b^J for transforming the position and velocity vectors in x' , y' , and z’ (see fig. 3) 

;o components in X,Y,Z is needed. This transformation matrix is identical with the 
ja^jJ transformation matrix where o)j, S2j, and ij are replaced by o>, S2, and i. 

The velocity vector at t = tj in the x' direction is p^x’), y' is p^(y'), z T 
is p i ( z') where 


P^x') = p a cos e i - p i e i sin 0^ 
%( y') = ^ sin 0 t + p i 8 i cos e i > 
Pi(z') =0 


(B42) 


The velocity components at t = tj in the x' ,y' ,z' coordinate directions are 


Pj(x') = p f cos 0 f - p f 0 f sin 9^ 
Pf(y’) = Pf sin 0 f + p f 0 f cos 9f > 
p f (z’) = 0 


(B43) 


These velocity components may be transformed to components in the X,Y,Z coordinate 
system by use of the (byl matrix above. 


From these results the velocities may be determined at 
ceptor orbit and t = tf in the target plane. Then at t = tj 
system: 


iiW = if eos "i - sin v i 


t = tj in the initial inter- 
in the x,y,z coordinate 


lj(y) = ^ sin v i + ZfVf cos v i 

li(z) = 0 


(B44) 


which may be transformed to X,Y,Z 


coordinates by the matrix 
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Finally, obtain the velocity after the final impulse directly from the components 
X,Y,Z without transformation 


• • m 

r^(X) = r f cos sin <p^ 

• • • V 

r f (Y) = r f sin + r^j cos <p^ / 
r f (Z) = 0 


(B45) 


The Velocities in X, Y, and Z 

The matrix operations on these vectors give the desired velocity components in the 
Newtonian frame X,Y,Z. These components are denoted by V k (j) where the subscript 
k refers to the time (where k equals 1 and 2 just before and after the initial impulse 
and k equals 3 and 4 just before and after the final impulse) and j to the component 
X,Y,Z. These velocity components, in terms of the components shown in equations (B42), 
(B43), (B44), (B45), and the elements of the [jayj and Q>ij] matrices, are given in the 
flow diagram of appendix D. 

The velocity increments required to perform the rendezvous maneuver are easily 
found by subtracting these components. For the velocity increment at t = t^, it is neces- 
sary to subtract the corresponding components of state (1) from state (2) as 


AV i (X) = V 2 (X) - V 1 (X) 

AV^Y) = V 2 (Y) - V 1 (Y) 

AVj(Z) = V 2 (Z) - V 1 (Z) 

and for the terminal maneuver subtract state (3) from state (4) to give the desired 
increments 


AV f (X) = v 4 (x) - v 3 (x) 
AV f (Y) = V 4 (Y) - V 3 (Y) 
AV f (Z) = V 4 (Z) - V 3 (Z) 


The total velocity increment for the maneuver is the sum of the maneuvers at t = t^ ^AVj 
and t = tj f AV 2 ^ 
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AV f = ^AV f (X ) 2 + AV f (Y ) 2 + AV f (Z ) 2 

and 

AV = AVj + AV f 
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DESCRIPTION OF PROGRAM 

This program was written in the FORTRAN IV (Ibsys version 9) language for the 
IBM 7094 computer at the Langley Research Center. Throughout the program there is an 
emphasis on simplicity, but capability has been provided for several uses and for a free- 
dom of choice on input and output. 

The method of solving the problem can readily be obtained by following the flow dia- 
gram (appendix D) and the description given in appendix B. The description in this 
appendix gives additional information about the options available, types of input and out- 
put, criteria for testing variables for transfer in the program, criteria for testing the 
convergence in the iterative processes, and criteria for incrementing the times. The 
methods used in the iterative processes are described fully in appendix B and can be 
followed on the flow diagram (appendix D). The flow diagram also shows the methods of 
incrementing At w , At, and v Q . 

A complete listing of the FORTRAN IV program is given as appendix E. 

Subprograms 

Six small subprograms are used in addition to the main program. They are called: 
AAA, VELOC, SIX, ANOM, SPACE, and CONVRT. An explanation of the computations, 
provided by AAA and ANOM, are included in appendix B. Equations (B5) are contained in 
AAA and equations (B9) are contained in ANOM. The other subprograms are self- 
explanatory; however, a brief description of the uses of all six subprograms is given here. 

AAA is used to compute the eccentric anomaly if the true anomaly is known or vice 
versa. 

VELOC is used to compute the velocity components. 

SIX is used to compute the elements of the Ja^ and [bj-y matrices. 

ANOM is used to compute the eccentric anomaly if the mean anomaly is known. 

SPACE is used to test the line count and to skip pages and print column headings when 
necessary. 

CONVRT is used to convert the times from dimensionless quantities to units corre- 
sponding to input for printing purposes. 


44 



APPENDIX C 


Options Available 

The program offers several options for choosing the variable which will be incre- 
mented initially and for choosing the type of input. This choice is made by reading in 
three control factors: GUIDE, OPTION, and DAT AS. 

GUIDE determines whether At w or At varies initially in the program and gives 
the appropriate output format regarding the choice. The program provides that the time 
At^y or At not incremented initially will be incremented secondly, after which v Q will 
be incremented. This procedure works for any number of At, At^, and v values. 

The quantity incremented initially will appear as the abscissa for ease in plotting or ana- 
lyzing the data. Hence, 

if GUIDE = 1, At w varies initially; 
if GUIDE = 2, At varies initially. 

OPTION provides a choice of three types of input for times: 


OPTION = 1, 

At w> 

At 

input dimensionless; 

OPTION = 2, 

At w> 

At 

input in minutes; and 

OPTION = 3, 

At w> 

At 

input in days. 


Restrictions are placed on maximums for At w and At because of the six spaces 
allowed by the output format for printout of the quantity varying initially. The restrictions 
are due only to the output format, and may easily be changed. 

(a) If either time exceeds 9999.0 minutes, they must be read in either dimensionless 
or in days. 

(b) If either exceeds 9999.0 days, they must be read in dimensionless. 

(c) If either exceeds 9.9999 in dimensionless time and is restricted by (a) or (b), the 
leftmost characters of the value will not be printed on output. 

DATAS provide a choice of input: 

if DATAS = 1, use Keplerian input; 
if DATAS = 2, use relative orbital input. 

Input 

Information on input may be found in the "Computer Input and Output" section of the 
paper and in the immediately preceding paragraphs. Input is to be made in units 
according to the following criteria: 
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(1) The quantities [i and arp are dimensioned, unless time is dimensionless in 
which case they may be in any units desired. However, in order to get values for P T 
and V^rp, ju. and a T should be dimensioned. 

(2) All angles are in degrees. 

(3) Times are in units described under "Options Available" in this appendix. 

(4) All other quantities are dimensionless. 

Input cards in the order to be read in and the proper FORTRAN formats for each are 
listed in the following table: 


Order of cards 

Variable names 

FORTRAN format 

1 

GUIDE, 0PTI0N, DAT AS 

4E18.8 

2 

AMU, AT 

4E18.8 

3 

AH, RMIN, ET, PHIO 

4E18.8 

4 

DTWI, DTWF, DTTI, DTTF 

4E18.8 

5 

DELTW, DELTT 

4E18.8 

6 

Keplerian input: 

El, AI, 0MEGAI, 0MI 

4E18.8 


Relative orbital input: 
XO, YO, ZO 

4E18.8 

7 

Keplerian input: 

AH, ANUOl, ANUOL, DELNUO 
Relative orbital input: 

XODT, YODT, ZODT 

4E18.8 


Output 

All output will be in units corresponding to input. The output format will vary as At 
and At w are varied initially. When a set of data is read in, those initial conditions are 
printed. As v Q and At or At w are incremented, they are printed. The computed 
values are then printed in columns. Special notation and messages which may be printed 
are 

(1) If relative orbital input is used and aj is computed as less than zero, a message 
to this effect will be written and the program will transfer to the initial input section. 

(2) If i = 90° - 0.1°, the message "i = +/-90 degrees not acceptable data” is 
written. 

(3) If the orbit is parabolic, the value for a will appear as 99.999 instead of °°. 
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(4) If p min < r min , an asterisk is printed at the extreme right-hand end of the line 
of data. 


Testing Criteria 

There are several places in the program where tolerance factors or allowances for 
computational inaccuracies must be defined to insure proper flow through the program. 
An explanation of the values chosen in these tests follows: 


Values compared 

Remarks 

h o3 

with 0.99999995 
H o 

This test compares the value of the cos ip to and 

restricts the minimum value of a computed ij to be 0.0001 
radian; otherwise, ij is set equal to 0 

1 i 1 with 90° - 0.1° 

Data which result in an i of 90° are not acceptable in this pro- 
gram. The absolute value of i is tested against 90° with a 
margin of 0.1° 

ATp with At 

This test of the parabolic time against the transfer time to 
determine the type of orbit results in a parabolic orbit only 
if the values are equal. A tolerance of 0.5 X 10“ 6 (At) defined 
as CRIT was allowed at the point of equivalence for computa- 
tional error 

i with 0.5 X 10- 7 

For this test a margin of tolerance of 0.5 X 10“ ^ was allowed 
for machine inaccuracy 

AT with At 

This test is made to determine when the convergence is suffi- 
cient to leave the loop in the iteration process. A value 
CRIT = 0.5 X 10"6(At) W as defined to give a tolerance mar- 
gin based on the value of the transfer time and to cover any 
computational error 

60 with 1.0 X 10" ^ 
6e with 1.0 x 1CT 7 

The appropriate test is made in each iteration scheme to test 
for the effectiveness of the increment on the variable. If the 
increment is equal to or less than the value tested against, 
its effect on the variable is negligible and the iteration is 
ended. This procedure acts as a safety check for ending the 
iteration in the event that AT never gets within the pre- 
scribed range (CRIT) of At 

Pj with p 2 

A tolerance margin of 1.0 X lO" 6 ^, defined as CRIRH0, pro- 
vides for transfer to the special iteration necessary when 
Pl=P 2 

A 0 with 0.0001° 
A 0 with 360° 

- 0.0001° 

Restrictions are put on a margin of 0.0001° around 

A 0 = 0°(360°) because of computational sensitivity. For 
these cases, transfer orbit properties are set equal to the 
initial target orbit properties and the iterations are omitted 
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The quantities used for testing in incrementing At w , At, and v Q are as follows: 
TWFMDT = At w (last) - 1.5(6tJ - TESTTW 
TTFMDT = At(last) - 1.5(6t) - TESTTT 
ANLMDN = t» 0 (last) - 1 . 5 ( 81 ^ - 1.0 X 10" 7 
DTWFT = At w (last) - TESTTW 
DTTFT = At(last) - TESTTT 


ANUOLT = j/ (last) - 1.0 X 10" 7 

TESTTW = 0.5 X 10 _6 jAt w (last)J, or = 1.0 X 10" 7 

TESTTT = 0.5 X 10" 6 [At(last)], or = 1.0 X 10" 7 


The TWFMDT is used to test against At w and allows it to be incremented by 6t, 


w 


until At w is greater than or equal to At w (last) - 1.5 6t w . A value TESTTW is arbitrar- 


w 


w* 


ily chosen as 1.0 x 10" 7 or 0.5 x 10“®[At w (last)| depending upon whether 6t w is 0 or is 
greater than 0, respectively. The TESTTW value is included to provide a tolerance for 
computational error in the value of TWFMDT if the 6t w is equal to 0 and to insure 
that control is transferred out of this loop after computation with At w (first). The DTWFT 

is used to test against 
At w (first) 


At. 


At w (last) 


w 


and always provides the computation for At w (last), even if 


6t 


w 


is not an integer. This separate control for the final point also 


insures ease in incorporating any additional desired statements or controls, such as those 
for plotting routines, in the program at the end of a series of incrementations. The 
TESTTW in the expression for DTWFT provides tolerance for machine error. 


The same explanation applies for TTFMDT, DTTFT, and TESTTT, if these variable 
names replace TWFMDT, DTWFT, and TESTTW in the preceding paragraph. These vari- 
ables are used for testing in the incrementing procedure for At. 

A similar explanation for ANLMDN and ANUOLT, used for incrementing v , holds 
with the exception that the tolerance margin or computational error is chosen as 
1.0 x 10" 7 for both values, instead of being a ratio as in TESTTW and TESTTT, because 
v is in radians (in the form X.XXXXXXX) and the 1.0 x 10" 7 will always be effective on 
the number of places used in the computation. 
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The mathematical symbol and its FORTRAN equivalent are given in the following table 


Mathematical symbol 

Fortran symbol 

a 

AT ,AI, A 

a ij 

A1 1 ,A1 2, A21 ,A22,A3 1 ,A3 2 

b ij 

B11,B12,B21,B22,B31,B32 

c 

C 

e 

ET,EI,E 

E 

ETO,ETF,EIO,En,EL,E2 

F 

F1,F2 

h 

AH 

H 

AH01,AH02,AH03,AH0 

i 

bi,ah,i,ii 

M 

AMTO,AMTF,AMIO,AMn,AMl 

P 

PI,PT,P,PP 

P T 

TPER 

r 

RO,RODT,RF,RDTF,RDTFl,RDTF2,RMIN 

At w 

DTWI,DTWF,DELTTW 

At 

DTTI,DTTF,DELTTT 

T 1> T 2 

T1,T2,TP1,TP2 

AT 

DT,DELTTP,DELTTC 

V 

V01,V02,V03,V0 

V CT 

VCT 

v k 0) 

V11,V12,V13,V21,V22,V23,V31,V32,V33,V41,V42,V43 

AV i (j),AV f (j) 

DV11,DV12,DV13,DV21,DV22,DV23 

AV i7 AV f 

DELV1.DELV2 

AV 

DELV 

x TT ,y' f ,z ,t ,x ,, ,y tf ,z” 

XO,YO,ZO,YODT,ZODT 

a 

ALPHA 

y 

GAMMA 

y 

GAMBAR 

A 6 

DELTTH 

6 

TH1,TH2,THP,THI,THF,THDTI,THDTF,THPMAX 

66 

DELTH 

P 

AMU, MU 

V 

ANUOl ,ANUOL, ANU0,ANUI,ANUDTI,NU01 ,NUOL,NUO,NUI 

k 

XIOl ,XI02,XI03 ,XIO,XIODT ,XH,XIDTI,Xni ,xn2,XH3,XnDT 1 , 
XHDT2 

P 

RH01,RH02,RH0I,RH0F,RH0DTI,RH0DTF,R0DTI1,R0DTI2, 

R0DTF1,R0DTF2,RH0M 

<P 

PHIO,PHIODT,PHIF,PHmTF 

* 

PSI 

00 

0MI,0M 

ft 

0MEGAI,0MEGA 
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in this diagram refers to Function Subprogram AAA. 
* + S2 in this diagram refers to Subroutine Subprogram ANOM. 
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Matrix M same as M 

with n replacing 

c i ) replacing coj 
i replacing ij 


i 


iiW = 

ij cos 

v i ' 

^ji^j sin ^ 

, 

i 

4i(y) = 

ij sin i^j + £ji/j cos i- 

i. 

i 

Pj(x’) = 

Pi 

COS 

9i 

" Pi^i 

sin 

h 

Pi(y') = 

Pi 

sin 

h 

+ pj0i 

cos 

8i 

Pf(x') = 

Pf 

COS 

6 f 

- Pf0f 

sin 

e f 

Pf(y') = 

Pf 

sin 


+ Pf0f 

cos 

Of 

H?* 

g 

11 

*f 

cos 

0f 

- r f 0f 

sin 

<£f 

r f (Y) = 

f f 

sin 

0f 

+ rf<£f 

cos 

<k 

nr - 








Vj(X) = anij(x) + a 12 ii(y) 

L 


V 3 (X) = bnp f (x’) +bi2Pf(y’) 

Vj(Y) = a 2 iii(x) + a 22 ij(y) 



V 3 (Y) = b2ipf(x’) +b 22 p f (y’) 

Vl(Z) = a 31 ij(x) + a 32 ij(y) 



V 3 (Z) = b 3 ip f (x’) + b 32 pf(y') 

V 2 (X) = bnpj(x’) + bi2Pi(y') 



V 4 (X) = r f (X) 

V 2 (Y) = b2lpj(x’) + b22Pi(y f ) 



V 4 (Y) = r f (Y) 

V 2 (Z) = b 3 ipj(x’) + b 3 2Pi(y') 



V 4 (Z) = 0 
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AVi(X) = V2(X) - Vi(X) 
AVi(Y) = V 2 (Y) - Vi(Y) 
AVi(Z) = V 2 (Z) - V^Z) 

AVf(X) = V4(X) - V 3 (X) 
AVf(Y) = V 4 (Y) - V 3 (Y) 
AVf(Z) = V 4 (Z) - V 3 (Z) 


I 


AVi = ^AVi(X) 2 + AVi(Y) 2 + AVi(Z) 2 
AVf = ^AVf(X) 2 + AVf(Y) 2 + AVf(Z) 2 
AV = AVi + AV f 


i 


Convert radians to degrees for printing: <fif, iq, A 6, i 

Convert times (At and At w ) into units corresponding to input 

1 


//rite computed quantities according to appropriate format 

j i : ~ ~ 


Convert times (At and At w ) back to dimensionless times 


I 
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FORTRAN PROGRAM 

The FORTRAN program to determine the velocity increment required for rendez- 
vous between two arbitrary elliptic orbits is as follows: 

WRITE(6,104) 

URITE(6,102) 

J = 7 
N J = 2 

10 READ ( 5, 100 ) GU I DE, OPT I OU, DATAS 
READ (5, 100) AMU, AT 
READ (5,100)/\!!, RJ 1 1 fl, ET,Pi:l 0 
READ ( 5, 100 )DTU I , DTWF, DTT I , DTTF 
READ (5,100) DELTW, DELTT 
I F(HJ.GT.1)G0 TO 48 
J = 3 

GO TO 49 

48 J = J 

49 !J J =0 

I F ( DATAS. OT. 1.5) GO TO 11 
R E A D ( 5 , 1 0 0 ) E I , A I , 0 ! 1 E G A 1 , 0 M I 
READ ( 5 , 1 00 ) A I I > AMUO 1 , AMU0 L, DELNU0 
24 T P E R = G . 2 8 3 1 8 5 3 * S QR T ( ( AT * * 2 / A.’ i 'J ) * AT ) 

vct=sqrt(a:iu/at) 

UR I TE ( G, 108 ) 

WRI TE(G, 106 )ET, El ,AI ,0.".EGAI ,0(11 , A! 1 
UR ITE(6,110) 

UR I TE( G, 106 ) Pii I 0 , A!!U0 1, AI.’UOL, Ai :U, AT, TPER, VCT 
UR I TE( 6, 112 ) 

UR I T E ( 6 , 1 2 6 ) DTT I , DTTF , RTW 1 , 0 TUF , R: 1 1 N , A! 1 
J = J+ 6 

I F (DATAS .GT. 1 . 5 )GO T025 
UR 1TE(6,124) 

J=J + 2 

55 UR I TE ( 6, 132 ) AIIU01 

C IF DTT I IS I i 1 PUT AS ZERO, IT IS REASS I OREL' A VALUE E<-d I VALENT TO 

c Ti :e ii!cre:;e!!t size. 

I F(DTT! .GT.0. 0)GO T r i 5 4 
IS IG = 4 9 
CALL SPACE 
J=J + 2 

DTT I = DELTT 
WRITE (6, 162) DTT I 
54 COIlTiniJE 

I F ( G'J I DE . GT . 1 . 5 ) GO TO 30 
I F (O PT I Oil . LT . 1 . 3 ) GO TC 50 
UR IT E( 6,140) DTT 1 
GO TC 51 

5 0 UR I TE( G, 130 ) DTT I 

51 UR I TE( 6, 116 ) 

GO TO 3 3 

30 COUTinUE 

I F (OPT I Of! . LT. 1 . 5 )GO TO 53 

UR I TE ( G, 133 ) ! jT!.'I 
GO TO 52 

53 UR I TE ( 6, 128 ) DTi; I 

52 UR I TE(6, 118 ) 

31 J=J+4 

I F (DATAS. GT. 1.5) GO TO 2G 
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PI=AI*(1.0-EI**2) 

PT“1. 0-ET**2 

26 PHi0=PHI0*.17453293E-01 

OMEGA I =OMEGAI *. 17453293E-01 
OMI =0MI *. 17453293E-01 
Ai I =A I I *. 17453293E-01 
ANU01=ANU01*. 17453293 E- 01 
ANUQ L=ANU0 L* . 17453293E-01 
DELNU0=DELNU0*. 17453293E-01 
GO TO 200 

11 READ (5/ 100 ) XO, YO, ZO 

READ(5/100)X0DT/Y0DT/Z0DT 
DELNU0=0 . 0 

PHI O-PH 10*. 174532 93 E-01 
PT=1 . 0-ET**2 
CPH I 0 = C0S ( PH 1 0 ) 

SPH I 0=S I N (PH I 0 ) 

R0=PT/(1.0+ET*CPHI 0) 

PH I ODT =SQRT (PT) / ( R0**2 ) 

RODT=(ET*SPHI 0)/SQRT(PT) 

V01=R0DT+X0DT-PH I ODT*YO 
V02=Y0DT+PH IODT*(RO+XO) 

V03=Z0DT 
X I 01 = R0+X0 
X I 0 2=Y0 
X I 03=Z0 

X I 0=SQRT(X IQ1**2+XI02**2+XI03**2) 

V0=SQRT(V01**2+V02**2+V03**2 ) 

AI=X10/(2.0-XI0*V0**2) 

C IF AI IS LESS THAN OR EQUAL TO ZERO/ A MESSAGE TO THIS EFFECT 1 

C WRITTEN AND TRANSFER IS TO THE BEGINNING OF THE PROGRAM/ 

C WHERE A NEW SET OF DATA MAY BE READ IN. 

1 F C AI >17/17/18 

17 TPER=S.2831353*SQRT( (AT**2/AMU ) *AT) 

VCT=SQRT( AMU/AT) 

WRITE (6/ 148) 

WR I TE ( 6/ 1 0 6 ) PH I 0 / ET/ AMU/ AT/ TPER/ VCT 
WRITE (6/ 112) 

WRITE(6/126)DTTI/DTTF/DTWI/DTWF/RMI U,R)j 
WR I TE( 6/ 114 ) 

WRITE(6/106)XO/YO/ZO/XODT/YGDT/ZODT 
VIR I TE( 6 / 124 ) 

WRITE (6/ 144) 

WR I TE (6, 104 ) 

GO TO 10 

18 AH01-XI02*V03-XI 03*V02 
AH02-X I03*V01-XI 01*V03 
AH03=X I 01*V02-X 1 02*V01 

PI =AH01**2+AH02**2+AH03**2 
AHO=SQRT(PI ) 

XI 0DT=(V01*XI01+V02*XI 02+V03*X I 03 )/X I 0 
PI DAI =PI /AI 
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TEST TO INSURE THAT El IS NOT EQUAL TO THE SQRT OF A NEC QUANTITY 
DUE TO SLIGHT COMPUTING INACCURACY ON PI/AI. ASSIGN VALUES 
FOR COS(NUO) AND SIN(NUO) IF THE QUANTITY UNDER THE RADICAL 
IS NEGATIVE. 

I F ( P I D A I . LT. 1. 0 )GO TO 19 
E I = 0 . 0 
CNU0=1 . 0 
SNU0 = 0 . 0 
GO TO 28 

19 CONTINUE 
El =SQRT( 1.0- (PI/AI ) ) 

CNU0= CCPI/X10) — 1.0)/EI 
SNUO=(SQRT(PI ) )*X1 ODT/E I 

28 CONTINUE 

ANU0=ATAN2 (SNUO, CNUO ) 

I FCSNUO.LT. 0.0. AND. CNUO. GT.O.O)AH'J0=6. 2831353 +ANU0 
I F (SNU 0 . LT. 0 . 0 . AND . Cf iU 0 . LT . 0 . 0 ) ANUO =6 . 2 8 3 1 8 5 3+ ANU Q 
13 CONTINUE 
ANU01=ANU0 
ANUO L=ANU0 
15 Cl I =Ali03/A!!0 

I F ( C 1 I ,GE. 0.99999995)00 TO 16 
A I I =ATAN (SQRT Cl. 0 — C I I **2)/Ci I ) 

I F ( C I I . LT . 0 . 0 ) A I I =3. 1415927+A! I 
SAI I =S I N ( A 1 I ) 

COMEG I =-(AH02*CPH I 0+AH01*SP! MO)/ ( ANO*SA 1 I ) 

SOM EG I = ( AHO 1*CPH I 0-AH02*S PH I 0 )/ ( AHO*SA I I) 

OMEGA ! =ATAN2 (SOMEG I , COMEG I ) 

I F (SGMEG I .LT. 0.0. AMD. COMEG I . GT. 0 . 0 ) OMEGA I =6 . 2831353+QMEGA! 

I F (SOMEG I . LT. 0.0. AND. COMEG I . LT. 0 . 0 )OMEGA I =6 . 28318 53+CMEGAI 

20 CONTINUE 

CO PN= ( X I 0 2 *AH01-X I 0 1*AH0 2 ) / ( X I 0*AH0 *SA I I ) 

SQPN=XI03/(XI0*SAI I ) 

GO TO 21 
IS A I 1=0.0 
OMEGA I =0. 0 

COPN=(XI 01*CPNI 0 — X I 02*SPHI0)/XI 0 
S0PN=(XI01*SPHI O+XI 0 2 * C PH I 0 ) / X I 0 

21 0PM=ATAN2 (SOPN, COPN) 

I F (SOPN. LT. 0. O.AND.COPN.GT. 0. 0)OPN=6. 2831853+OPN 
I F (SOPN. LT. 0.0. AND. COPN. LT. 0.0 )0PN=6. 28318 53+0 PM 
23 CONTINUE 

ONI I =OPN-ANUO 

C CONVERT ANGLES IN RADIANS TO DEGREES FOR PRINTING 

PH I 0 = PHI 0/ . 1745329 3E-01 
OMEGA I =OMEGA I / . 17453293E-01 
OM I =014 I/. 174 53293 E- 01 
Al I =AI I/.17453293E-01 
ANU01=ANU01/ . 17453293E-01 
ANUOL=ANUOL/ .1745 3293 E-01 
DELNU0=DELNU0/ . 17453293 E- 01 
GO TO 24 
25 WR I TE ( 6/ 114 ) 

WRITE(6,106)XO,YO,ZO,XODT,YODT,ZODT 
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WRITE (6/ 12 4) 

J=J + 4 
GO TO 55 

C CHANGE TIMES TO DIMENSIONLESS QUANTITIES 

200 I F(QPTI ON. LT. 1 . 5 ) GO TO 201 
I F (OPT I ON. LT .2. 5) GO TO 202 
GO TO 203 

202 CM I N=60 . /TPER 
DTW I =DTW I *CM I N 
D T W F = D T W F * CM I N 

D E LTW = D E LTW * C M I N 
DTTI =DTTI *CMI N 
DTTF=DTTF*CM I N 
DELTT=DELTT *CM I N 
GO TO 201 

203 CDAYS=86400 ,/TPER 
DTWI=DTWI*CDAYS 
DTWF =DTWF *CDAYS 
DELTW=DELTW*CDAYS 
DTTI =DTTI *CDAYS 
DTTF=DTTF*CDAYS 
DELTT=DELTT*CDAYS 

201 CONTINUE 

I F (DTWF. GT. 0.0)60 TO 204 
TESTTW*. 0000001 
GO TO 205 

204 TESTTW*. 5E-06*DTWF 

205 I F (DTTF . GT. 0 . 0 )GO TO 206 
TESTTT b . 0000001 

GO TO 207 

206 TESTTT=.5E-06*DTTF 

207 CONTINUE 

TWFMDT =DTWF-1 . 5*DE LTW-TESTTW 
TTFMDT-DTTF-1 . 5*DE LTT-TESTTT 
AN LMDN=ANU0L~1. 5 *DELNU0- .0000001 
DTW F T= DTW F - T ES TTW 
DTTFT=DTTF-TESTTT 
ANU0LT-ANU0L-. 0000001 
I F(GU I DE-1. 5)39,43/43 

39 ANU0=AMU01 

40 DELTTT=DTTI 

41 DELTTW=DTW I 
GO T042 

43 ANU0=ANU01 

44 DELTTW=DTWI 

45 DELTTT=DTTI 

C CRIT IS A VALUE USED TO TEST THE ITERATION SCHEMES FOR CONVERGENCE 

42 CR I T= . 5E-0 6*DE LTTT 
CR I RHO=l . 0 E-06*A I 
COMMONJ^GUI DE / ISIG 
GAMMA" 1 . 0 

61 ET0=AAA(ET, PHI 0, GAMMA) 

AMT0=ET0-ET* (S I N(ET0 ) ) 

AMTF =AMT0+6 .2831853* (DELTTT+DE LTTW ) 
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CALL ANOM ( ET, AliTF , ETF ) 

63 PK I F=AAA(-ET , ETF , GAMMA ) 

CPK I F =COS ( PH I F ) 

SPHI F=SI N(PHI F) 

RF=PT/(1.0+ET*CPHI F) 

PHI DTF=(SQRT(PT) )/(RF**2) 

RDTF=ET*SPH I F/SQRT ( PT) 

65 E I 0 =AAA (El, ANUO, GAMMA ) 

AM I 0 = E I O-E I * (S I N ( E 1 0 ) ) 

AMI I =AM I 0+ ( ( 6 . 2831853*DE LTTW) / (Al**1.5)) 

CALL ANOM ( E I , AM I I , E 1 1 ) 

67 ANU I =AAA(-E I , E I I, GAMMA) 

XI i=PI/(1.0+EI *COS ( ANU I )) 

ANUDT I =SQRT ( P I )/ ( X 1 I **2 ) 

X I DT I =E I *S I N (AMU I ) / (SQRT ( P I )) 

C DETERMINATION OF DELTTH AND I OF TRANSFER ORBIT 

300 0 I PN I =ANU I +OM I 

CO I PNI =COS(OI PNI ) 

SOI PNI -SI N(OI PNI ) 

CO MEG I =COS(QMEGAI ) 

SOMEGI -SIM (OMEGA I ) 

Cl 1 =COS (A1 l ) 

SI l-SINCAl I ) 

XI 1 1=X I I * ( COMEG I *CO I PN I - C I I *SOMEG I *SO I PN I ) 

XI 1 2 =X I I * (SOMEG I *CO I PN I + C I I *COtlEG I *SC I PN I ) 

XI I3=XI I* (SI I *SO I PNI ) 

SUB1=X I I 3/ (X I I 1*SPII I F-X I I 2*CPH IF) 

B I =ATAN(SUB1 ) 

ASOLI =ABS ( B I ) 

I F(ASOLI -1.5690510 >306,305, 305 

305 CONTINUE 
I S I G=48 
CALL SPACE 
WR 1 TE ( 6, 150 ) 

J= J+2 

GO TO 1028 

306 CONTINUE 

CDTH= ( X I 1 1*CPH I F+X I I 2*SP!i I F )/X I I 

SDTH= (X I 1 1*S PH I F-X I I 2 * C PH I F ) / ( X I l*COS(BI )) 

DELTTH=ATAN2 (SDTH,CDTH) 

I F (SDTH. LT . 0 . 0 . AND . CDTH . GT . 0 . 0 ) DELTTH = 6 . 28 31853+DE LTTI i 
IF(SDTH. LT. 0.0. AND.CDTH.LT. 0.0) DELTTH=6. 28 31853+DELTTli 
302 CONTINUE 

C COMPUTATION OF PARABOLIC TIME 
RM02 =AMAXl ( RF , X I I ) 

RH01=AM I N1 ( RF, X I I ) 

I F (DELTTH. GT.O. 1745E-04.AND.DELTTH.lt. 6. 2831679)GO TO 308 
DEGDTH=DELTTH/ , 17453293 E- 01 
I S I G=47 
CALL SPACE 
J=J + 1 

WRI TE(6, 160 )DEGDTH 

A-1.0 

P= PT 
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TH1=PH I F-DELTTH 
TH2=PH I F 
E=ET 

GO TO 800 
308 CONTINUE 

C=SQRT (RH01**2+RH02**2-2 . 0*RH01*RH02*CDTH) 

SPS I = (Rr!01-RH02*CDTH)/C 
C PS I = ( RHO 2 *S DTH ) / C 
PSI=ATAN2(SPSI / CPSI) 

IF (SPS I .LT. 0.0. AND. CPS I . LT. 0 . 0 ) PS I =6. 2831853+PS I 
304 CONTINUE 

SALP=(RH02-RH01 )/C 
I F(SALP.GE.1.0)G0 TO 311 
ALPHA=ATAN(SALP/SQRT (1. Q-SALP**2 ) ) 

GO TO 312 

311 ALPHA=1. 5707963 

312 CONTINUE 
THP-ALPHA-PSI 
PP=RHGl* ( 1. Q.+ COS (THP) ) 

TTH P=(S I NCTHP/2 . ) ) / ( COS (T; IP/2 . ) ) 

TTHPDT=(S I i!( (THP+DELTTIO/2. ))/(OS((TuP+DELTTH)/2. )) 

TP1=( (PP**1.5)*( (TTHP**3)/6.+TTI!P/2 . ) )/ 6. 28 318 53 
TP2=( (PP**1 . 5 ) *( (TTHPDT**3 )/6 , +TTNPDT/2 . ))/G. 2851853 
DELTTP=TP2-TP1 

THPMAX=3 . 14159 2 7- (ALPHA* PS I ) 

SU312=DELTTP-DELTTT 
ASUB12=ABS (SUB12 ) 

I F ( ASUB12-CR I T) 37/ 37/ 3 8 

C IF DELTTP EQUALS DELTTT,THE ORBIT IS PARABOLIC 

37 E=1 . 0 
P=PP 
TH1=THP 

TH 2 = TH 1 + D E LTTH 

C A IS INFINITE/ IS SET EQUAL TO 99.939 FOR PRINTING PURPOSES 
A=99. 999 
GO TO 800 

C IF DELTTP IS GREATER THAN DELTTT/TNE ORBIT IS HYPERBOLIC 

C IF DELTTP IS LESS THAN QELTTT/THE ORBIT IS ELLIPTIC 

38 I F (DELTTP-DE LTTT) 400/ 37/ 600 

400 SUB13=RHQ1-RHQ2 
ASUB13=ABS (SUB13 ) 

I F (ASUB13-CR 1 R-HO >500 / 500/ 401 
C ELLIPTIC LOOP RH02 NOT EQUAL TO RHOl 

401 DELTH=(3. 1415927-2. *(THP+PSI ))/AH 
TH1=THP+DELTH 

402 TH2=TH1+DELTTH 

j^E= ( RH02-RH01 ) / (RH01*C0S (TH1 ) - R! 10 2* COS (TH2 ) ) 

^Di FF-THPMAX-TK1 

I F ( D I FF . LE . 0 . 0 ) GO TO 405 
I F ( E . GE . 1 . 0 )GO TO 405 
P=RH01* (1 . 0+E*COS (TH1 ) ) 

A=P/(1.0-E**2) 

GAMMA=1 .0 

E1=AAA(E/TH1/ GAMMA ) 



APPENDIX E 


E 2 = AAA ( E / TH 2 , G Ai li 1A ) , 

DT=((A**1.5)*(E2-El-E*(SirKE2)-SIN(El))))/6.2831853 

SUB14=DT-DELTTT 

ASUB14=ABS (SUB14 ) 

IF(ASUB14-CRIT)800,300/403 

403 IF(SUB 14 >404/800/405 

404 TH1=TH1+DELTH 
GO TO 402 

405 DELTH=DELTH/2. 

I FIDELTH-. 0000001)800/ 800/ 408 
408 TH1=TH1-DELTH 
GO TO 402 

ELLIPTIC LOOP RH02 EQUAL TO RHOl 

500 DELTTC=(DELTTH*(RH01**1. 5 ))/6. 2831853 
DELE=1.0/AH 

E=DELE 

SUB15=DELTTC-DELTTT 
IF(SUB 15)502/501/501 

501 P=RH01* (1.0+E*(C0S(PSI ))) 

A=P/ (1 . 0-E**2 ) 

TH1 — 1.0*PSI 
TH2=TH1+DELTTH 
GAMMA = 1 . 0 

E 1 =A AA ( E / TH 1 / GAMMA ) 

DT=- C(A**1.5) * (E1-E*S I NCEl) ) )/3. 1415927 

SUB16=DT-DELTTT 

ASUB16=ABS (SUB16 ) 

I F(ASUB16-CR!T)800/800/503 

503 IF(SUB16)504/800/505 

504 DELE=DELE/2. 

I F(DELE-. 0000001)800/ 800, 512 
512 E=E-DELE 
GO TO 501 

505 E=E+DELE 
GO TO 501 

502 P=RHOl* (1 . 0-E* (COS ( PS I ))) 

A=P/(1.0-E**2) 

TH 1 = 3. 1415927- PS I 
TH2=TH1+DELTTH 
GAMMA =1.0 

E1=AAA(E/TH1/ GAMMA) 

AMl=El-E*(SliK(El) ) 

DT=( (A**l. 5 ) *(3. 1415927-AMI) )/3. 1415927 

SUB17=DT-DELTTT 

ASUB17=ABS(SUB17) 

IF(ASUB17-CRIT)800,800/506 

506 I F(SUB17)507/800,508 

507 E=E+DELE 
GO TO 502 

508 DELE=DELE/2 . 

I F(DELE-. 0000001)800, 800, 514 
514 E=E-DELE 
GO TO 502 
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600 SUB18=RH01-RH02 
ASUB18=ABS(SUB18 ) 

I FCASUB18-CRI RHO 1700,700/ 601 
C HYPERBOLIC LOOP RH02 NOT EQUAL TO RHOI 

601 DELTH=(THP+PS I )/AH 
610 TH1=THP-DELTH 

60 2i^TH2 =TH1+D E LTTH 

"r e = (RH02-RH01)/( RHOI* (COS (TH1))-RH02* (COS (TH2))) 
CTH1=C0S (TH1 ) 

OPECT=1.0+E*CTH1 
I F(OPECT.GT.O.O)GO TO 608 
609 DELTH=DELTH/2. 

I F (DELTH . GT.0.1E-08)GO TO 612 


C 

C 

C 


611 

612 

608 


IF THE CASE IS A LIMITING HYPF.i-.BOL I 
(-1.0 )/COS (DELTTH/2 . 0 ) '..'RITE A 
COMPUTATION FOR THAT CASE. 

I S I G=49 
CALL SPACE 


ORBIT l 
MESSAGE 


HTi! E APP 
AND DISCO: 


J = J + 2 

E- (-1. 0 )/COS (DELTTH/2. 0) 

I F (GU I DE . GT. 1 . 5 ) GO TO 611 

CALL COIIVRT(DELTTM / CDTl.' / OPTICn / TPER) 

WRITE(6, 164)CDTW, E 
GO TO 1028 

CALL COIIVf;T(DELTTT,CDTT, OPTI ON, TPER) 

WR I TE( 6, 166 ) CDTT, E 
GO TO 1028 
CONTINUE 
TH1=TH1+DELTH 
GO TO 602 

I F(E . LE . 1 . 0 )GO TO 604 
P=RHC1*(1.0+E*(CGS(TN1)) ) 

A=P/ ( 1 . 0-E**2 ) 

GAM; 1A— 1.0 
F 1 = AA A ( E , TH 1 , G A. IMA) 

F2=AAA(E,T!!2,GAiiMA) 

TANF 1=(SIN (?!))/( COS (FI)) 

5U319=Fl/2. +3. 1415027/4. 

TS19 = (Siri(SUB19))/(CCS(SUG13)) 

Tl = ( ((-A)**l. 5) *(E*TAKF1-AL03(TS12 )) )/C. 2831853 
TANF2 = (SIN(F2))/( COS ( F2 ) ) 

SUS20=F2/2. +3. 1415027/4. 

TS20= (S I N (SU02C ; )/ (CCS (SJ020 ) ) 

T2= ( ( (-A) **1 . 5 ) * (E*TA!iF2 -ALCG (TS2 0 } ) )/ 6 . 28 318 53 

DT-T2-T1 

SUB21=DT-DELTTT 

AS U C 2 .1 = AO S ( S UP 21) 

I F ( ASUB21- CR IT) 8 00 ,8 00 ,5 03 
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•30 3 

1 F(SU:;21) 

GO 4 

:)Elt:>:>el' 


i foelt::- 

:o7 

Cfi.TIfi'JE 
t; ,l=y' : i p 


-***•> ~ r> 

J k t* '.j — 

305 

70 GO?. 



700 

J2LE = 1 ,/.V 


c ;.l to 


:=i.+ ."Li 

y;;i =- ps | 

CT;:i=cuScn:i ) 

t!.2=t::i-:-deltt;; 

701 p=n!:oi*(i.+E*C 7 ::i) 

A=P/ (1. ~E**2 ) 

.0 

i : 1=A/.A (E ,7\ \i, 

TF1=(S I i i ( F 1 ) ) / ( COG (FI) ) 

SUL522-F1/2. + 3. 14 15927/4. 

TS22*(£ I :iCS'Ji;22 ) )/ (HOG (GU'22) ) 

T =-(((- A ) * * 1 . 5 ) * ( E * 7“l-ALv C7S22)))/ 
S U 3 2 3 = PT - 2 E L TTT 
A SUE 2 3= AES (SUF. 23 ) 

I F ( AS U I ; 2 3 - 0 : : 1 T ) C 0 0 , 8 4 0 , 7 0 2 

702 I F(SU 123)703,8 00,70 4 

703 PE LE = EELF./2 . 

IF (DELE-. 0000001)800,800, 70 7 
707 E-E-JELE 
GO TO 701 

704 E=E+USLE 
GO TO 701 

000 I F( TADS 01,802,802 

coi a::c;;k;=p/(i.+e) 

30 TO 003 
00 2 P.Ii-.i ; I L = ;V :01 
GC T,: 30 3 
203 :c= 0 

I F ( KI :c. i!;:-A!H f:)S04,C05, 005 
304 K=2 

805 I F ( ilF - X I I >306/307/ 3 0 7 
006 R;:Oi =R!iG2 
ii::cF=, r u:ci 

T::i =6. 2831853-T;:2 
TKF =Ti ! 1 -i'LELTT! : 

GO TO 303 
30 7 :>HC, I = 11:01 

r;:of=ri;g2 

T!i I =Ti!l 
T!!F=T!';2 

308 $QP=SQP,T(P) 

THDTI=Sf t P/(P.!!OI**2) 

P.KOOT 1 =E*S I N (TIM )/SQP 
Ti!OTF=SQP/(RI!OF**2) 
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. RHODTF =E*S I N (THF ) /SC-P 

/ /VI x~rhl2 ,-t&T7P2Z7K5T, A3 2 , 01 11,0: IEGA I ,A1D 

809 A3I-ABSCQI ) 

IF(AB 1-5. £-08)810,810, 811 

810 OMEGA-O. 0 


on=pi:iF-Tf;F 
GO TO 900 

811 I F ( 3 I )8 13, 810,812 

812 0 * 1 EGA " PIIIF+3. 1415927 


813 

900 


Oi'i=3. ltl59 27~T!'F 
GO TO 900 
Oi1EGA=Pli I F 


0: i=5 . 28 318 53-THF 


CALL SIX(C 11, 012,821,8 22, 33 1,032,0.' 1, 0! iEGA, B I ) 

XI I DT1=X ! DT 1 * ( COS ( AMU I ) )-Xi 1 *Af!UDTI *(S I IKAMUI ) ) 
X I I DT2 =X I DT i * (S 1 :i(Al!U I ) ) + X I 1 *ANUDT I * ( COS (AN U I ) ) 


RODT1 l=RHODTI * (COS (TIM ) )-R!!C I *THDTI *(SI IKTHI )) 
ROOT! 2=RHGQTI *(S I N(TM I ) ) + RHC ! *Tf!DTI * (COS (Till ) ) 
RODTFl = Rl!ODTF*(COS(T!IF) )-RI!QF*THDTF*(S I H(THF) ) 
RODTF2 = RHODTF* (S I N (THF ) ) + RIIOF *THDTF * ( CCS (TKF ) ) 


R0TF1«RDTF*( COS (Pill F) )-RF* PH I DTF* (S I fj (P1M F) ) 

RDTF2 =RDTF* (S I M ( PH I F ) ) + RF * PH I DTF* ( COS ( PH 1 F ) ) 

901 CALL VE LOG (Vll, V12, V13, All, A12, A2 1, A2 2, A31, A32 , X I : DTI, XI IDT?) 
CALL VELOC(V21,V2 2, V2 3, C11,B12,B21, 322,531, B32, P.ODT II, ROOTI 2) 
CALL VELQC(V31, V32, V33, Bll, 812, 821, B22, D31, B32, RGDTF 1, R0DTF2 ) 
V41=RDTF1 


V42=RDTF2 
V43=0 . 0 

1001 DV11=V21-V11 
DV12=V22-V12 
DV13=V2 3-V13 
DV21=V41-V31 


DV22=V42-V32 


DV23=V43-V33 

1002 0ELV1=S0RT(DV11**2 + DV12 **2 + DV13**2 ) 
DELV2=SQRT(DV21**2+DV22**2+DV23**2) 
DELV=DELV1+DELV2 
PHI F=PH IF/. 174532 93E-01 
AND I =ANU I / . 174532S3E-01 
TH I =TH I / . 174532 93 E-0,1 
DELTTH-DELTTH/. 174532 93E-01 
BI=BI/.17453293E-01 

CALL CGNVRT ( DE LTTT, DE LTTT, OPT I OH, TPER) 
CALL C 0 !J V R T ( D E L TTl , D E L TTW , 0 P T i 0 N , T P E R ) 
I S I G = 49 
CALL SPACE 
J=J +2 


I F(GUIDE.GT.1.5)GO TO 1030 
IF(K.GT.l)GO TO 1032 
I F (OPTI ON . GT. 1 . 5 ) GO TO 1040 
WR !TE(6,120)DE LTT'W, PI; I F, AMU I , E , A, 3 
113, DV21,DV22, DV23,DELV1,DELV2,DELV 
GO TO 1031 






72 



APPENDIX E 


1040 URITE(6,134)DELTTW.,PH|F,ANUI,E,A,3I ,1m I ,DELTT!I, RMC. H N, DV11, DV12, OV 
113,DV21,DV22,DV23,DELV1,DELV2, DELV 

GO TO 1031 

1032 I F (OPT I ON. GT. 1.5) GO TO 1041 

WRI TE( 6/ 12 2 ) DELTTW, PH I F , AMU I , E, A, B 1 , TIM , DE LTTII, RIIG; 1 1 N, DV1 1, DV12 , 3V 
113,DV21,DV22,DV23,DELV1,DELV2,DELV 
GO TO 1031 

1041 UR I TEC 6,136 )DELTTU, PH | F,ANU.I , E, A, G 1 , TH I , DELTT!:, PNOM I N, DV11 , DV1 1’ DV 
113,DV21,DV22,DV23,QELV1,DELV2,DELV 

GO TO 1031 

1030 I F(K.GT.l)GO TO 1033 

I F(OPTI QN.GT . 1 . 5)GO TO 1042 

UR I TE( 6 / 12 0 )DE LTTT , Pit I F, AMU I , E, A, 8 1 , Til I , DELTTH, RIIO! i ! N, DV11, DV12 , DV 
113,DV21,DV22,DV23,DELV1,DELV2,DELV 
GO TO 1031 

1042 WRI TE (6/134) DELTTT , Pi ! I F, ANU I , E, A, G I , TNI ,QELTT!!, RNQi 11 N, DV11, DV12, DV 
113, DV2 1, DV2 2,DV2 3, DE LV1, DE LV2, DELV 

GO TO 1031 


1033 1 F (OPT lOfl. GT. 1. 5 )GO TO 1043 

VJR I TE ( 6 , 1 2 2 ) D E LTTT, PI i I F / ANU I / E / A, B ! / T! ! I / D E LTTi ! 
113/DV21/ DV2 2/ DV2 3/ DELV1, DELV2/ DE LV 
GO TO 1031 

104 3 UUITE(6/13G) DELTTT / PI 1 1 F/AIHJ I , E , A, D I , T! : I , DELTT!! 

113/ DV21/0V22/ DV23/DELV 1/DE LV 2/DE LV 
1031 COflT I HUE 


nnp ' ; 1 H 
/ t \ 1 ; V 1 I I < / 

/ iihOi ! I i' 1 / 


3 VI 1, DV12 , DV 


C CONVERT SACK TO DIMENSIONLESS TIME 

IF(GPTI0M.LT.1.5)G0 TO 1020 
I F (CPTI Ofl. LT. 2 . 5 )GQ TO 1072 
DELTTT=3G400 . *DE LTTT/TPER 
DELTTU=3 6400 , *DELTTU/TPER 
GO TO 1028 


1072 ,DELTTT = 60. *DE LTTT/TPER 
DELTTU=GO.-*DELTTU/TPER 
1023 CONTINUE 


IF (GUI DE-1. 5) 1003/ 1015/ 1015 

1003 I F(GELTT'.:-TUF:DT)1004/1C05/ 1005 

1004 0 E LT Tl .' = D E L TTU + D E LTI 
N0 = 0 

GO TC42 

1005 I F (NO . GT. 1)GG TO 1007 1 ' 

I F ( D E LTT! L, T : F T ) 1 J06/ 1 007,1007 

1006 DELTTU=DTUF 
NO = 2 

GO TO 42 

1007 I F(DE LTTT-TTF; IDT) 1003 / 1000 / 100 0 
1003 DE LTTT=R E LTTT *DE LTT 

* 10=0 


1101=0 
I S I G=48 


CALL SPACE 
URI TE(6/ 124) 

CALL CQNVRT ( D E LTTT, CDTT, 0 PT 1 ON, TP E ' : ) 
I F(CPTI ON . GT. 1 . 5 )GO TO 10G0 
WA! TE( 6/130 )C!;TT 
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GO TO 1061 

1060 WRITE(6,140)CDTT 

1061 J=J + 3 

GO TO 41 

1009 IF(NQ1.GT.1)G0 TO 1011 

I F(DELTTT-DTTFT) 1010 / 1011/ 1011 

1010 DELTTT =DTTF 
N 0 = 0 

NO 1=2 
ISI G=48 
CALL SPACE 
WRI TE ( 6/ 124 ) 

CALL CONVRT(DELTTT/CDTT/OPTl ON/TPER) 
IF(OPTIOH.GT.1.5)GO TO 1062 
WRITE(6,130) CDTT 
GO TO 1063 

1062 WRI TE (6, 140 )CDTT 

1063 J=J+3 

GO TO 41 

1011 I F (ANUO-ANLMDN )1012,1013/1013 

1012 AM U 0 =AN U 0 + D E LHU 0 
N0=0 

N01=0 
N02 = 0 
1 S I G = 4 7 
CALL SPACE 
WRITE(6/124) 

ANU0=AMUp/ .17453293E-01 
WR I TE( 6/ 132 ) AMUO 

CALL CONVRT(DTT I / CDTT I /OPT I ON/TPER) 
IF(OPTION.GT.1.5)GO TO 1064 
WRI TE (6/ 130 ) CDTT I 
GO TO 1065 

1064 WR I TE ( 6, 140 )CDTT I 

1065 ANU0=ANU0* . 174532 93 E-Ol 
J=J + 3 

GO TO 40 

1013 I F(N02.GT.1)G0 TO 12 

I F(ANUO-ANUOLT) 1014/ 12/12 

1014 ANU0=ANU0L 
N0=0 

NO 1 = 0 
NO 2 =2 
I S I G=4 7 
CALL SPACE 
WR I TE ( 6/ 124 ) 

ANUO-ANUO/ . 17453293E-01 
WRITE (6, 132) AMUO 

CALL CONVRT ( DTT I / CDTT I /OPTION/TPEP.) 

I F (OPT ION. GT. 1,5) GO TO 1066 
WRITEC6/130) CDTT I 
GO TO 1067 

1066 WRI TE ( 6 , 140 ) CDTT I 

1067 ANU0=ANU0*. 17453293E-01 
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J=J + 3 
GO TO 40 

1015 I F (DELTTT-TTFi IDT) 1016 / 1017/ 1017 

1016 DELTTT=DELTTT+DELTT 
NO-O 

GO TO 42 

1017 I F (NO.GT. DGQ TO 1019 

I F(DELTTT-DTTFT) 1018/ 1019/1019 

1018 DELTTT=DTTF 
NO = 2 

GO TO 42 

1019 I F (DELTTW-TWF! IDT) 1020 / 1021/ 1021 

1020 DELTTW-DELTTW+DELTW 
H0=0 

noi-o 

I S I G=48 
CALL SPACE 
WR I TE ( 6/ 12 4 ) 

CALL CONVRTCDELTTW, CDTVi/OPTI ON/ TPER) 
IF(OPTIOI1.GT.1.5)GO TO 1080 
WRI TE ( 6/ 12 8 ) CDTW 
GO TO 1081 

1080 WRI TE( 6/ 138 ) CDTW 

1081 J =J + 2 

GO TO 45 

1021 IF(N01.GT.1)G0 TO 1023 

I F (DELTTW-DTWFT) 1022, 1023/ 1023 

1022 DELTTW=DTUF 
NO=0 
N01=2 

I S I G=48 
CALL SPACE 
WRI TE( 5, 124 ) 

CALL COflVRTfDELTTV.'/ CDTW/ 0 PT I ON/ TPER) 
I F (OPTI ON. GT. 1 . 5 >GO TO 1082 
WRITE(6,128 ) CDTW 
GO TO 1083 

1082 V/R I TE ( 6 / 13 8 ) CDT'ij 

1083 J=J + 2 

GO TO 45 

1023 I F(ANUO-ANLMDN) 1024/ 1025/ 1025 

1024 ANUO=ANUG+DELNUO 
NO = 0 

NO 1=0 
NO 2=0 
I S I G=47 
CALL SPACE 
WR I TE ( 6/ 12 4 ) 

ANU0=ANU0/. 174532 93E-01 
W R I TE ( 6 / 1 3 2 ) ANU 0 

CALL CONVRTCDTWI /CDTWI / OPT I ON/ TPER ) 

I F (OPT ION. GT. 1.5) GO TO 1084 
WRITEC 6/128 )CDTW 1 
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GO TO 1085 

1084 WRI TE ( 6, 138 ) CDTW I 

1085 ANU0=ANU0*. 17453293E-01 
J=J+3 

GO TO 44 

1025 I F (N02 . GT. 1 )G0 TO 12 
IF!ANU0-ANU0LT)1026,12,12 

1026 ANU0=ANU0L 
NO=0 
N01=0 

NO 2 =2 
I S I G=47 
CALL SPACE 
WRITE! 6, 124) 

ANUO-ANUO/. 174532 93E-01 
WRI TE( 6, 132 ) ANUO 

C A L L C 0 N V R T ( D T W I , C D T W I , 0 P T I 0 N , T P E R ) 

I F (OPT I ON .GT . 1. 5 )GO TO 1086 
WRLTE(6,128) CDTW I 
GO TO 1087 

1086 WR I T E ( 6 , 1 3 8 ) CDTW I 

1087 ANU0=ANU0*. 17453293E-01 
J=J + 3 

GO TO 44 
12 WRITE(6,104) 

GO TO 10 

100 FORMAT! 4 E18 . 8 ) 

102 F ORMAT ( / / 4 X 4 2 H V E LO C I T Y INCREMENT REQUIRED FOR RENDEZVOUS IX 3 7HSETWE 
1EN TWO ARBITRARY ELLIPTIC ORBITS//) 

104 FORMAT (1H1) 

106 FORMAT! 7E18. 8) 

108 FORMAT! 12X2HET16X2HE I 16X2HAI 12X6HOMEGAI 15X3HOMI 16X2NI I ) 

110 FORMAT! 10X4HPH I 014X4HNU0114X4HNU0 L16X2HMU16X2HAT14X4HTPER15X3HVCT) 
112 FORMAT! 10X4HDTT I 14X4HDTTF14X4HDTW I 14X4HDTWF32X4HRHI N17X1HI!) 

114 FORMAT! 12X2HX0 16X2 HY016X2HZ0 14X4HX0 DTI 4 X4HY0DT14X4HZ0DT) 

116 FORMAT! IX 6HDELTTW3X4H PH I F4X3NNUI 6X1HE6X1HA6X1H I 5X3HTH 1 1XGHDELTTH1X 
14HRHOM4X4HDV1X4X4HDV1Y4X4HDV1Z4X4HDV2X4X4HDV2Y4X4HDV2Z2X5HDELV12X 
25HDELV23X4HDELV//) 

118 FORMAT! 1X6HDE LTTT3X4H PH I F4X3HNU I 6X1HE6X1HA6X1H I 5X3HTH 1 1X6HDELTTH1X 
14HRHOM4X4ilDVlX4X4HDVlY4X4HDVlZ4X4HDV2X4X4HDV2Y4X4HDV2Z2X5HDE LV12X 
25HDELV23X4HDELV// ) 

120 FORMAT! F 7 .4,2F7.2,F7.4,F7.3,F7.2,F8.2,F7.2,F5.2,6F8.5,3F7.4//) 

122 FORMAT !F7.4,2F7.2, F7.4,F7.3,F7.2,F8.2,F7.2,F5.2,6F8.5,3F7.4,1H*//) 
124 FORMAT ! IX// ) 

126 FORMAT !4E18.8,E36.8,E18.8) 

128 FORMAT! IX 8HDELTTW= F7.4) 

130 FORMAT ! 1X8HDE LTTT - F7.4) 

132 F ORMAT ! 4X5HNUO= F7.2) 

134 FORMAT! F7. 1,2F7.2,F7.4,F7.3,F7.2,F8.2,F7.2,F5.2,6F8.5,3F7.4//) 

136 FORMAT! F 7 . 1,2F7.2,F7.4, F7.3,F7.2,F8.2,F7.2,F5.2,6F8.5,3F7.':,1H*//) 
-138- -FORMAT -C 1-X&4D E LT-TW *= F7 T l-> 

140 FORMAT! 1X8HDELTTT= F7.1) 

144 FORMAT! 6X20HA I IS LESS THAN ZERO//) 

148 FORMAT ! 10X4HPH I 016X2HET16X2HMU16X2HAT14X4HTPER15X3HVCT) 
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150 FORMAT(4X41HI EQUAL +/-90 DEGREES NOT ACCEPTABLE DATA//) 

160 FORMAT (10X19HDELTA THETA EQUALS F6. 1,1X39 H DEGREES. DO NOT CONS I DE 
1R UNLESS ORB I TS1X25H INTERSECT, IN WHICH CASE,) 

162 FORMAT C 6X4 8 HDTT I INPUT EITHER AS ZERO OR A NEGATIVE QUANT I TY1X23HB 
1UT EQUATED TO DELTT ( F10.4,1X33H) BY PROGRAM. RECONSIDER INPUT. 
2 //) 

164 FORMAT ( 10X20HTH I S CASE (DELTTW = F10.4,1X37H) IS A LIMITING HYPERB 
10LIC ORBIT W I TH1X14HE APPROACHING F10.5//) 

166 FORMAT ( 10X20HTH I S CASE (DELTTT = F10.4,1X37H) IS A LIMITING NY PERB 
10LIC ORBIT W I TH1X14HE APPROACHING F10.5//) 

END 


FUNCTION AAA C AL PH A , B ETA , GAI IMA ) 

GAMMA EQUAL 1, EXCEPT FOR HYPERBOLIC ORBITS WHEN IT IS -1 
P I NUM=3. 1415927 
BETA2=BETA/2 . 0 
ABETA2=ABS ( BETA2 ) 

SBETA2=S I N(BETA2 ) 

I F(SBETA2.GT. 0.0. AND. ABETA2.GT.1. 5707789. AND. ABETA2.LT.1~. 570 8137 )G 
10 TO 1 

I F(SBETA2. LT. 0.0. AND. ABETA2.GT. 1.5707789. AND. ABETA2.LT. 1. 5708137 )C- 
10 TO 2 

SUB 1= (SQRT (GAMMA* ( 1 . 0 -ALPHA)/ ( 1 . 0+ ALPHA) ))*(SBETA2/(COS(BETA2)) ) 
APR=2 . *ATAN (SUB1 ) 

GO TO 3 

1 APR=3. 1415927 
GO TO 3 

2 APR=-3. 1415927 

3 N=(BETA+PI MUM)/ ( 2 . *PI HUM) 

AN=N 

AAA=APR+2 . *AN*PI HUM 

RETURN 

END 


SUBROUTINE VELOCC Yll, Y12, Y13, H11,H12, H21, H22 , H31, H32 ,S1, S2 ) 

Y11=H11*S1+H12*S2 

Y12=H2 1*S-1 «U22a$2 

Y13=H31*S1+H32*S2 

RETURN 

END 
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SUBROUT I HE S I X(W11,W12/W21 / U22 / 1731/W32/ U1/U2/ U3 > 
CU1=C0S(U1) 

SU1-SIMCU1) 

CU2=COS ( U2 ) 

SU2=S I M(U2 ) 

CU3=COS (U3 ) 

SU3=S I N(U3 ) 

WU=CU1*CU2-CU3*SU2*SU1 

U12 =-SUl*CU2-CU3*SU2*CUl 

W2 1 = CU 1 *S U2 + CU 3 * CU 2 *S U 1 

W22=-SU1*SU2+CU3*CU2*CU1 

W31=SU3*SU1 

W32=SU3*CU1 

RETURN 

END 


SUBROUTINE ANOM(ECCEN, AMANOM, EANOH) 

BARM-AMANOM 

E-ECCEN 

E AO =BA RM+E*S I fl (BARM )+.5*E**2*SIN(2. *B ARM ) 
EA=EAO 

16 AMA-EA-E*SIH(EA) 

DELMA-BARM-AMA 

DE LEA=DELMA/ (1. -E*COS (EA) ) 

EA1=EA+DELEA 

ADELEA=ABS(DELEA) 

CONTST® (EA1-EA)/EA1 
ACT =ABS ( CONTST) 

1 F (ACT- .1E-06)17/15/15 
15 EA=EA1 
GO TO 16 

17 EANOM-EA 
RETURN 
END 


SUBROUTINE SPACE 
COMMON J, GU I DE, I S I G 
IF(ISIG-J)10/10/11 
10 WR I TE( 6/ 101 ) 
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I F (GU I DE . GT . 1 . 5 )G0 TO 12 
WR I TE( 5, 10 3 ) 

GO TO 13 

12 WRITE(6 / 104) 

13 J = 5 


11 

101 

103 


104 


RETURN 
CONTINUE 
RETURN 
FORMAT QH1) 


FORMAT ( 1X5HDELTTW3X4HPI i I F4X3NNU I 6X1NE6X1KA5X1N I 5X3NTH i 1X5IDE LTTK1X 
14HRHOM4X4HDV1X4X4HDV1Y 4X4HDV1Z4X4HDV2X4X4HDV2Y4X4HDV2Z2X5HDELV12X 
25HDELV23X4NDELV//) 

FORMAT ( 1X6ND ELTTT3X4NP1-I I F4X3MNU I 6X1NE6X1HA0X1N ! 5X3HTN i IX 6 HOE LTTH1X 
14HRH0M4X4HDVlX4X4HDVlY4X4HDV124X4i;0V2X4X4!i0V2Y4X4!DV2Z2X5!:0ELY12X 


25HDELV23X4NDELV// ) 
END 


SUBROUTINE C0NVRT(T I ME, CT I ME, OPT I ON, TPER) 

TO CONVERT TIME FROM DIMENSIONLESS QUANTITIES TO UNITS 
CORRESPONDING TO INPUT 
IFCOPTIOM.LT. 1.5)GO TO 1 
IFCOPTION.LT. 2. 5)G0 TO 2 
CTI ME=T IME*TPER/8 6400 . 

RETURN 

1 CT I i 1 E =T I ME 
RETURN 

2 CT I ME=T IME*TPER/ 60 . 

RETURN 

END 
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TABLE I 

COMPUTER PRINTOUT OF EXAMPLE 1 


VELOCITY INCREMENT REQUIRED FOR RENDEZVOUS BETWEEN TWO ARBITRARY ELLIPTIC ORBITS 


ET 

El 

A I 

OMEGA I 


UMI 

II 



0. 933 71999E-01 

0.16724200E-01 

0 . 65630 IOOE 00 

0.25387999E 

03 

0 .23 3C1999E 03 

0* 18499999E 01 



PHIO 

NUO I 

NUOL 

MU 


AT 

TPER 

VCT 


0 • 324 39999E 03 

0. 37C000CCE-00 

0.3 70CC0QQ E-00 

0. 46789999E 

22 

0.7 4 736999E 12 

0.59348101E OB 

0.79124083E 

05 

DTT I 

DTTF 

DTW I 

DTWF 



RMIN 

H 


C. 16000CCCE 03 

0.2600C000E 03 

0 . 

0.400CGOOOE 

02 


0.39378100E-00 

0. 30000000E 

01 


NUO- 0.37 


DELTTW = 
DELTTT 

= C, 

PHI F 

NUI 

E 

A 

. 

THI 

DELTTH 

RHOM 

DV IX D VI Y 

DV 1Z DV2X DV2Y 

DV2Z 

DELV1 

DELV2 

OELV 

160.0 

423.41 

C .37 

0.6758 

0.627 

1.65 

-134.33 

296.14 

0.20 

0.78032-0.44667 

0.03413-0.51167-0.05421 

0.01771 

0.8998 

0.5148 

1.4146* 

18C.0 

434.75 

0.37 

0.6661 

0.596 

1.87 

-136.86 

307.48 

0.20 

0.76728-0.40968 

0.02693-0.44464 0.01959 

0.01944 

0.8702 

0.4455 

1.3157* 

200.0 

445.71 

0.37 

0.6383 

0.583 

2.24 

-139.61 

318.43 

0.18 

0.80293-0.40036 

O.G1994-0. 44264 0.10123 

0.02198 

0.8974 

0.4546 

1.3520* 

220. C 

456.29 

0.37 

0.7463 

0.581 

2.38 

-143.67 

329. CO 

C . 15 

0.90147-0.40946 

0.01285-0.51051 0.17498 

0.02549 

0.9902 

0.5403 

1.5305* 

240. C 

466.52 

0.37 

C . 8392 

0.586 

4.19 

-150.66 

339.22 

0.09 

1.C8507-0. 41653 

0.00574-0.65745 0.21178 

0.02985 

1.1623 

0.6914 

1.8537* 

260.0 

476.43 

0.37 

0.9427 

0.598 

7.86 

-162.05 

349.07 

0.03 

1.36525-0.37671- 

-0.00035-0. 87692 0.15690 

0.03412 

1.4163 

0.8915 

2.3078* 

DELTTH 5 

160.0 

- 2C. 

434.75 

,C 

20.74 

0.6640 

0.668 

1.86 

-129.46 

287.10 

0.22 

0.87931-0. 192b2 

0.02853-0.46880-0.20692 

CD 

N»- 

o 

CM 

O 

O 

0.9006 

0.5128 

1.4135* 

180.0 

445.71 

2C.74 

0.6432 

0.630 

2.01 

-131.76 

298.06 

0.22 

C. 84048-0. 16181 

0.02048-0.40512-0.12722 

0.02168 

0.8562 

0.4252 

1.2813* 

200.0 

456.29 

20.74 

G. 6486 

0.612 

2.27 

-134.02 

308.63 

0.22 

C. 84514-0. 14257 

0 .01297-0. 39640-0.05278 

0.02360 

0.8572 

0.4006 

1.2578* 

220.0 

466.52 

20.74 

0.6817 

0.606 

2.70 

-136.89 

318 .86 

0.19 

0.89663-0.12717 

0.C0561-0. 43889 0.00647 

0.02635 

0.9056 

0.4397 

1.3453* 

240. C 

476.43 

20.74 

0.7448 

C *607 

3.43 

-141.31 

328.75 

0.16 

1.00350-0.10339-0.00183-0.53243 0.03607 

o 

o 

u> 

o 

o 

1.0088 

0.5345 

1.5433* 

26C.0 

486.05 

20.74 

0.8348 

0.614 

4.82 

-148.51 

338.34 

0.10 

1.17529-0.04870- 

•0.C0926-0. 67511 0.01106 

0.03454 

1.1763 

0.6761 

1.8524* 

DELTTW= 

160.0 

4C . 

445.71 

C 

41.03 

0.6567 

0.716 

1.86 

-124.49 

277.77 

0.25 

0.89817 0.07952 

0.02058-0.39029-0. 34103 

0.02106 

0.9019 

0.5187 

1.4206* 


DELTTT 

PHIF 

NUI 

E 

A 

I 

THI 

DELTTH 

RHOM 

DV IX 

DV1Y 0V1Z DV2X DV2Y 

DV2Z 

DELV1 

0ELV2 

DELV 

180.0 

456.29 

41, .03 

0.6258 

0.668 

1.94 

-126.59 

288.35 

C . 25 

0.83996 

0.09298 0 .01247-0.33246-0.25209 

0.02159 

0.8452 

0.4178 

1.2630* 

20C.0 

466.52 

41 .03 

C.6178 

0.645 

2.10 

-128.52 

298.57 

0.25 

0.81985 

0.10684 0.00515-0.32083-0.17773 

0.02273 

0.8268 

0.3675 

1.1943* 

220.0 

476.43 

41. C3 

0.6321 

0.634 

2.36 

-130.68 

308.47 

0.23 

0.83549 

0.126GO-0. 001 78-0. 34789-0. 12324 
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COMPUTER PRINTOUT OF EXAMPLE 2 
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COMPUTER PRINTOUT OF EXAMPLE 3 
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"The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof .” 

— National Aeronautics and Space Act of 1958 
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